Smart thermostat with model predictive control

ABSTRACT

A thermostat for a building zone includes at least one of a model predictive controller and an equipment controller. The model predictive controller is configured to obtain a cost function that accounts for a cost of operating HVAC equipment during each of a plurality of time steps, use a predictive model to predict a temperature of the building zone during each of the plurality of time steps, and generate temperature setpoints for the building zone for each of the plurality of time steps by optimizing the cost function subject to a constraint on the predicted temperature. The equipment controller is configured to receive the temperature setpoints generated by the model predictive controller and drive the temperature of the building zone toward the temperature setpoints during each of the plurality of time steps by operating the HVAC equipment to provide heating or cooling to the building zone.

CROSS-REFERENCE TO RELATED PATENT APPLICATIONS

This application is a continuation of U.S. patent application Ser. No.15/625,830 filed Jun. 16, 2017, which claims the benefit of and priorityto U.S. Provisional Patent Application No. 62/491,545 filed Apr. 28,2017. Both of these patent applications are incorporated by referenceherein in their entireties.

BACKGROUND

The present disclosure relates generally to a smart thermostat and moreparticularly to a smart thermostat with model predictive control.Thermostats are often configured to monitor and control the temperatureof a building zone or other space. For example, a thermostat can bemounted on a wall within the building zone and configured to measure thetemperature of the building zone. Thermostats typically send commands toHVAC equipment (e.g., on/off commands, heating/cooling commands, etc.)to cause the HVAC equipment to affect the temperature of the buildingzone.

Conventional thermostats operate according to a fixed temperaturesetpoint schedule which defines the temperature setpoints for thethermostat at various times. The temperature setpoint schedule istypically set by a user via a local user interface on the thermostat. Inmany implementations, a fixed temperature setpoint schedule leads tosuboptimal control of the HVAC equipment, which can increase the cost ofheating/cooling the building zone. It would be desirable toautomatically determine optimal temperature setpoints for a thermostatin order to take advantage of time-varying energy prices, zone heattransfer characteristics, and/or other factors that can affect the costof heating/cooling the building zone.

SUMMARY

One implementation of the present disclosure is a thermostat formonitoring and controlling temperature of a building zone. Thethermostat includes an equipment controller and a model predictivecontroller. The equipment controller is configured to drive thetemperature of the building zone to an optimal temperature setpoint byoperating HVAC equipment to provide heating or cooling to the buildingzone. The model predictive controller is configured to determine theoptimal temperature setpoint by generating a cost function that accountsfor a cost operating the HVAC equipment during each of a plurality oftime steps in an optimization period, using a predictive model topredict the temperature of the building zone during each of theplurality of time steps, and optimizing the cost function subject to aconstraint on the predicted temperature of the building zone todetermine optimal temperature setpoints for each of the plurality oftime steps.

In some embodiments, the model predictive controller is configured todetermine the cost of operating the HVAC equipment during each of theplurality of time steps using a set of time-varying utility ratescomprising a utility rate value for each time step. The time-varyingutility rates may be received from a utility provider or predicted bythe model predictive controller.

In some embodiments, the model predictive controller is configured topredict the temperature of the building zone during each of theplurality of time steps as a function of a temperature setpointtrajectory comprising a temperature setpoint for each of the pluralityof time steps.

In some embodiments, the model predictive controller is configured tooptimize the cost function subject to a constraint on the optimaltemperature setpoints that limits the optimal temperature setpointswithin a temperature setpoint range.

In some embodiments, the model predictive controller is configured togenerate the predictive model by performing a system identificationprocess. The system identification process may include modulating thetemperature setpoint within a constrained temperature setpoint range,collecting a set of input-output data, and fitting parameters of thepredictive model to the set of input-output data. The input-output datamay include values of the temperature setpoint and values of thetemperature of the building zone that result from modulating thetemperature setpoint during each of a plurality of time steps during alearning period.

In some embodiments, the predictive model includes a thermal massstorage model that defines the temperature of the building zone as afunction of at least one of heat transfer between air within thebuilding zone and solid mass within the building zone, heat transferbetween the building zone and the HVAC equipment, and an unmeasured heatload disturbance. In some embodiments, the model predictive controlleris configured to predict a value of the unmeasured heat load disturbanceexperienced by the building zone at each of the plurality of time stepsin the optimization period.

In some embodiments, the predictive model includes an HVAC load modelthat defines the heating or cooling provided by the HVAC equipment as afunction of the temperature of the building zone and the temperaturesetpoint.

In some embodiments, the model predictive controller is configured topredict the cost of operating the HVAC equipment as a function of anamount of the heating or cooling provided by the HVAC equipment at eachtime step of the optimization period.

In some embodiments, the constraint on the predicted temperature of thebuilding zone requires the model predictive controller to maintain thepredicted temperature of the building zone within a first zonetemperature range during a first time step of the optimization periodand within a second zone temperature range, different from the firstzone temperature range, during another time step of the optimizationperiod subsequent to the first time step.

Another implementation of the present disclosure is a method performedby a thermostat for a building zone for monitoring and controllingtemperature of the building zone. The method includes generating a costfunction that accounts for a cost operating HVAC equipment during eachof a plurality of time steps in an optimization period, using apredictive model to predict the temperature of the building zone duringeach of the plurality of time steps, optimizing the cost functionsubject to a constraint on the predicted temperature of the buildingzone to determine optimal temperature setpoints for each of theplurality of time steps, and operating HVAC equipment to provide heatingor cooling to the building zone to drive the temperature of the buildingzone to the optimal temperature setpoints.

In some embodiments, the method includes receiving a set of time-varyingutility rates from a utility provider or predicting the time-varyingutility rates. The set of time-varying utility rates may include autility rate value for each time step. The method may includedetermining the cost of operating the HVAC equipment during each of theplurality of time steps using the set of time-varying utility rates.

In some embodiments, using the predictive model to predict thetemperature of the building zone includes predicting the temperature ofthe building zone during each of the plurality of time steps as afunction of a temperature setpoint trajectory that includes atemperature setpoint for each of the plurality of time steps.

In some embodiments, optimizing the cost function includes optimizingthe cost function subject to a constraint on the optimal temperaturesetpoints that limits the optimal temperature setpoints within atemperature setpoint range.

In some embodiments, the method includes generating the predictive modelby performing a system identification process. The system identificationprocess may include modulating the temperature setpoint within aconstrained temperature setpoint range, collecting a set of input-outputdata, and fitting parameters of the predictive model to the set ofinput-output data. The input-output data include values of thetemperature setpoint and values of the temperature of the building zonethat result from modulating the temperature setpoint during each of aplurality of time steps during a learning period.

In some embodiments, the predictive model includes a thermal massstorage model that defines the temperature of the building zone as afunction of at least one of heat transfer between air within thebuilding zone and solid mass within the building zone, heat transferbetween the building zone and the HVAC equipment, and an unmeasured heatload disturbance. In some embodiments, the method includes predicting avalue of the unmeasured heat load disturbance experienced by thebuilding zone at each of the plurality of time steps in the optimizationperiod.

In some embodiments, the predictive model includes an HVAC load modelthat defines the heating or cooling provided by the HVAC equipment as afunction of the temperature of the building zone and the temperaturesetpoint.

In some embodiments, the method includes predicting the cost ofoperating the HVAC equipment as a function of an amount of the heatingor cooling provided by the HVAC equipment at each time step of theoptimization period.

Another implementation of the present disclosure is a thermostat formonitoring and controlling temperature of a building zone. Thethermostat includes an equipment controller and a model predictivecontroller. The equipment controller is configured to drive thetemperature of the building zone to a zone temperature setpoint byoperating HVAC equipment to provide heating or cooling to the buildingzone. The model predictive controller is configured to determine thezone temperature setpoint by generating a cost function that accountsfor a cost operating the HVAC equipment during each of a plurality oftime steps in an optimization period, using a predictive model topredict the temperature of the building zone during each of theplurality of time steps, and performing an optimization of the costfunction subject to a constraint on the predicted temperature of thebuilding zone to determine a temperature setpoint trajectory including atemperature setpoint value for each of the plurality of time steps.

Another implementation of the present disclosure is a model predictivecontrol system for monitoring and controlling temperature of a buildingzone. The model predictive control system includes a thermostat and amodel predictive controller. The thermostat is configured to drive thetemperature of the building zone to an optimal temperature setpoint byoperating HVAC equipment to provide heating or cooling to the buildingzone. The model predictive controller is configured to determine theoptimal temperature setpoint and provide the optimal temperaturesetpoint to the thermostat via a communications network. The modelpredictive controller determines the optimal temperature setpoint bygenerating a cost function that accounts for a cost operating the HVACequipment during each of a plurality of time steps in an optimizationperiod, using a predictive model to predict the temperature of thebuilding zone during each of the plurality of time steps, and optimizingthe cost function subject to a constraint on the predicted temperatureof the building zone to determine optimal temperature setpoints for eachof the plurality of time steps.

In some embodiments, the model predictive controller is configured todetermine the cost of operating the HVAC equipment during each of theplurality of time steps using a set of time-varying utility rates thatincludes a utility rate value for each time step. The time-varyingutility rates may be received from a utility provider or predicted bythe model predictive controller.

In some embodiments, the model predictive controller is configured topredict the temperature of the building zone during each of theplurality of time steps as a function of a temperature setpointtrajectory comprising a temperature setpoint for each of the pluralityof time steps.

In some embodiments, the model predictive controller is configured tooptimize the cost function subject to a constraint on the optimaltemperature setpoints that limits the optimal temperature setpointswithin a temperature setpoint range.

In some embodiments, the model predictive controller is configured togenerate the predictive model by performing a system identificationprocess. The system identification process may include modulating thetemperature setpoint within a constrained temperature setpoint range,collecting a set of input-output data, and fitting parameters of thepredictive model to the set of input-output data. The input-output datamay include values of the temperature setpoint and values of thetemperature of the building zone that result from modulating thetemperature setpoint during each of a plurality of time steps during alearning period.

In some embodiments, the predictive model includes a thermal massstorage model that defines the temperature of the building zone as afunction of at least one of heat transfer between air within thebuilding zone and solid mass within the building zone, heat transferbetween the building zone and the HVAC equipment, and an unmeasured heatload disturbance. In some embodiments, the model predictive controlleris configured to predict a value of the unmeasured heat load disturbanceexperienced by the building zone at each of the plurality of time stepsin the optimization period.

In some embodiments, the predictive model includes an HVAC load modelthat defines the heating or cooling provided by the HVAC equipment as afunction of the temperature of the building zone and the temperaturesetpoint.

In some embodiments, the model predictive controller is configured topredict the cost of operating the HVAC equipment as a function of anamount of the heating or cooling provided by the HVAC equipment at eachtime step of the optimization period.

In some embodiments, the constraint on the predicted temperature of thebuilding zone requires the model predictive controller to maintain thepredicted temperature of the building zone within a first zonetemperature range during a first time step of the optimization periodand within a second zone temperature range, different from the firstzone temperature range, during another time step of the optimizationperiod subsequent to the first time step.

Another implementation of the present disclosure is a method formonitoring and controlling temperature of a building zone. The methodincludes generating a cost function that accounts for a cost operatingHVAC equipment during each of a plurality of time steps in anoptimization period, using a predictive model to predict the temperatureof the building zone during each of the plurality of time steps,optimizing the cost function at a model predictive controller subject toa constraint on the predicted temperature of the building zone todetermine optimal temperature setpoints for each of the plurality oftime steps, and providing an optimal temperature setpoint from the modelpredictive controller to a thermostat for the building zone via acommunications network. The method further includes, at the thermostat,driving the temperature of the building zone to the optimal temperaturesetpoint by operating HVAC equipment to provide heating or cooling tothe building zone.

In some embodiments, the method includes receiving a set of time-varyingutility rates from a utility provider or predicting the time-varyingutility rates. The set of time-varying utility rates may include autility rate value for each time step. The method may includedetermining the cost of operating the HVAC equipment during each of theplurality of time steps using the set of time-varying utility rates.

In some embodiments, using the predictive model to predict thetemperature of the building zone includes predicting the temperature ofthe building zone during each of the plurality of time steps as afunction of a temperature setpoint trajectory that includes atemperature setpoint for each of the plurality of time steps.

In some embodiments, optimizing the cost function includes optimizingthe cost function subject to a constraint on the optimal temperaturesetpoints that limits the optimal temperature setpoints within atemperature setpoint range.

In some embodiments, the method includes generating the predictive modelby performing a system identification process. The system identificationprocess may include modulating the temperature setpoint within aconstrained temperature setpoint range, collecting a set of input-outputdata, and fitting parameters of the predictive model to the set ofinput-output data. The input-output data include values of thetemperature setpoint and values of the temperature of the building zonethat result from modulating the temperature setpoint during each of aplurality of time steps during a learning period.

In some embodiments, the predictive model includes a thermal massstorage model that defines the temperature of the building zone as afunction of at least one of heat transfer between air within thebuilding zone and solid mass within the building zone, heat transferbetween the building zone and the HVAC equipment, and an unmeasured heatload disturbance. In some embodiments, the method includes predicting avalue of the unmeasured heat load disturbance experienced by thebuilding zone at each of the plurality of time steps in the optimizationperiod.

In some embodiments, the predictive model includes an HVAC load modelthat defines the heating or cooling provided by the HVAC equipment as afunction of the temperature of the building zone and the temperaturesetpoint.

In some embodiments, the method includes predicting the cost ofoperating the HVAC equipment as a function of an amount of the heatingor cooling provided by the HVAC equipment at each time step of theoptimization period.

Another implementation of the present disclosure is a model predictivecontrol system for monitoring and controlling temperature of a buildingzone. The model predictive control system includes a thermostat and amodel predictive controller. The thermostat is configured to drive thetemperature of the building zone to a zone temperature setpoint byoperating HVAC equipment to provide heating or cooling to the buildingzone. The model predictive controller is configured to determine thezone temperature setpoint and provide the zone temperature setpoint tothe thermostat via a communications network. The model predictivecontroller determines the zone temperature setpoint by generating a costfunction that accounts for a cost operating the HVAC equipment duringeach of a plurality of time steps in an optimization period, using apredictive model to predict the temperature of the building zone duringeach of the plurality of time steps, and performing an optimization ofcost function subject to a constraint on the predicted temperature ofthe building zone to determine the temperature setpoint trajectory.

Another implementation of the present disclosure is a model predictivecontroller for monitoring and controlling temperature of a buildingzone. The model predictive controller includes a system identifier and apredictive optimizer. The system identifier is configured to operate themodel predictive controller in a system identification mode. Operatingin the system identification mode includes performing a systemidentification process to automatically generate a predictive modelbased on heat transfer characteristics of the building zone. Thepredictive optimizer is configured to operate the model predictivecontroller in an operational mode. Operating in the operational modeincludes using the predictive model to predict the temperature of thebuilding zone. The model predictive controller is configured toautomatically transition from the operational mode to the systemidentification mode in response to a determination that a predictionerror of the predictive model exceeds a threshold value.

In some embodiments, operating in the operational mode includesoptimizing a cost function subject to a constraint on the predictedtemperature of the building zone to determine optimal temperaturesetpoints for each of a plurality of time steps in an optimizationperiod. In some embodiments, the cost function accounts for a cost ofoperating HVAC equipment to provide heating or cooling to the buildingzone during each of the plurality of time steps in the optimizationperiod.

In some embodiments, performing the system identification processincludes collecting a set of input-output data and fitting parameters ofthe predictive model to the set of input-output data. The input-outputdata may include values of a temperature setpoint for the building zoneand the temperature of the building zone during each of a plurality oftime steps during a learning period.

In some embodiments, the set of input-output data includes a discreteHVAC staging trajectory for staged HVAC equipment. The HVAC stagingtrajectory may include a discrete HVAC equipment load at each of theplurality of time steps during the learning period. In some embodiments,performing the system identification process includes filtering theinput-output data to generate a continuous HVAC equipment load signalfrom the discrete HVAC staging trajectory.

In some embodiments, collecting the set of input-output data includesmodulating a temperature setpoint within a constrained temperaturesetpoint range and recording values of the modulated temperaturesetpoint and values of the temperature of the building zone that resultfrom modulating the temperature setpoint. In some embodiments,performing the system identification process comprises filtering theinput-output data to remove oscillations in the temperature of thebuilding zone around the temperature setpoint resulting from activatingand deactivating staged HVAC equipment.

In some embodiments, the predictive model includes a thermal massstorage model that defines the temperature of the building zone as afunction of heat transfer between air within the building zone and solidmass within the building zone.

In some embodiments, the model predictive controller includes aload/rate predictor configured to predict a value of an unmeasured heatload disturbance experienced by the building zone at each of a pluralityof time steps in an optimization period. In some embodiments, thepredictive model defines the temperature of the building zone as afunction of the unmeasured heat load disturbance.

In some embodiments, the predictive model includes an HVAC load modelthat defines an amount of heating or cooling provided by HVAC equipmentcontrolled by the model predictive controller as a function of thetemperature of the building zone and a temperature setpoint for thebuilding zone.

Another implementation of the present disclosure is a method formonitoring and controlling temperature of a building zone. The methodincludes operating a model predictive controller in a systemidentification mode. Operating in the system identification modeincludes performing a system identification process to automaticallygenerate a predictive model based on heat transfer characteristics ofthe building zone. The method further includes operating the modelpredictive controller in an operational mode. Operating in theoperational mode includes using the predictive model to predict thetemperature of the building zone. The method further includesautomatically transitioning from the operational mode to the systemidentification mode in response to a determination that a predictionerror of the predictive model exceeds a threshold value.

In some embodiments, operating in the operational mode includesoptimizing a cost function subject to a constraint on the predictedtemperature of the building zone to determine optimal temperaturesetpoints for each of a plurality of time steps in an optimizationperiod. In some embodiments, the cost function accounts for a cost ofoperating HVAC equipment to provide heating or cooling to the buildingzone during each of the plurality of time steps in the optimizationperiod.

In some embodiments, performing the system identification processincludes collecting a set of input-output data and fitting parameters ofthe predictive model to the set of input-output data. The input-outputdata may include values of a temperature setpoint for the building zoneand the temperature of the building zone during each of a plurality oftime steps during a learning period.

In some embodiments, the set of input-output data includes a discreteHVAC staging trajectory for staged HVAC equipment. The HVAC stagingtrajectory may include a discrete HVAC equipment load at each of theplurality of time steps during the learning period. In some embodiments,performing the system identification process includes filtering theinput-output data to generate a continuous HVAC equipment load signalfrom the discrete HVAC staging trajectory.

In some embodiments, collecting the set of input-output data comprisesmodulating a temperature setpoint within a constrained temperaturesetpoint range and recording values of the modulated temperaturesetpoint and values of the temperature of the building zone that resultfrom modulating the temperature setpoint. In some embodiments,performing the system identification process includes filtering theinput-output data to remove oscillations in the temperature of thebuilding zone around the temperature setpoint resulting from activatingand deactivating staged HVAC equipment.

In some embodiments, the predictive model includes a thermal massstorage model that defines the temperature of the building zone as afunction of heat transfer between air within the building zone and solidmass within the building zone.

In some embodiments, the method includes predicting a value of anunmeasured heat load disturbance experienced by the building zone ateach of a plurality of time steps in an optimization period. In someembodiments, the predictive model defines the temperature of thebuilding zone as a function of the unmeasured heat load disturbance.

In some embodiments, the predictive model includes an HVAC load modelthat defines an amount of heating or cooling provided by HVAC equipmentcontrolled by the model predictive controller as a function of thetemperature of the building zone and a temperature setpoint for thebuilding zone.

Another implementation of the present disclosure is a model predictivecontroller for monitoring and controlling temperature of a buildingzone. The model predictive controller includes a state/disturbanceestimator and a predictive optimizer. The state/disturbance estimator isconfigured to estimate an initial state of the building zone at abeginning of an optimization period. The state of the building zoneincludes a temperature of air within the building zone and a temperatureof solid mass within the building zone. The predictive optimizer isconfigured to use a predictive model of the building zone to predict thestate of the building zone at each of a plurality of time steps of theoptimization period based on the estimated initial state of the buildingzone and a temperature setpoint trajectory comprising a temperaturesetpoint for each of the plurality of time steps. The predictiveoptimizer is configured to generate optimal temperature setpoints foreach of the plurality of time steps by optimizing a cost function thataccounts for a cost operating HVAC equipment during each of theplurality of time steps.

In some embodiments, the predictive optimizer is configured to generatethe cost function and determine the cost operating the HVAC equipmentduring each time step using a set of time-varying utility rates thatincludes a utility rate value for each time step. The set oftime-varying utility rates can be received from a utility provider orpredicted by the model predictive controller.

In some embodiments, the predicted state of the building zone includes apredicted temperature of the building zone. The predictive optimizer maybe configured to optimize the cost function subject to a constraint onthe predicted temperature of the building zone.

In some embodiments, the model predictive controller includes aload/rate predictor configured to predict a value of an unmeasured heatload disturbance experienced by the building zone at each of theplurality of time steps in the optimization period. In some embodiments,the predictive model defines the temperature of the building zone as afunction of the unmeasured heat load disturbance experienced by thebuilding zone at each time step.

In some embodiments, the predictive model includes a mass storage modelthat defines the temperature of the building zone as a function of heattransfer between air within the building zone and solid mass within thebuilding zone. In some embodiments, the predictive model defines thetemperature of the building zone as a function of heat transfer betweenthe building zone and the HVAC equipment. In some embodiments, thepredictive model includes an HVAC load model that defines an amount ofheating or cooling provided by the HVAC equipment as a function of thetemperature of the building zone and the temperature setpoint.

In some embodiments, the predictive optimizer is configured to predictthe cost of operating the HVAC equipment as a function of an amount ofheating or cooling provided by the HVAC equipment at each time step ofthe optimization period.

In some embodiments, the model predictive controller includes a systemidentifier configured to generate the predictive model by performing asystem identification process. The system identification process mayinclude collecting a set of input-output data and fitting parameters ofthe predictive model to the set of input-output data. The input-outputdata may include values of the temperature setpoint and the temperatureof the building zone during each of a plurality of time steps during alearning period.

In some embodiments, performing system identification process includesmodulating the temperature setpoint within a constrained temperaturesetpoint range and recording values of the modulated temperaturesetpoint and values of the temperature of the building zone that resultfrom modulating the temperature setpoint.

Another implementation of the present disclosure is a method formonitoring and controlling temperature of a building zone. The methodincludes estimating an initial state of the building zone at a beginningof an optimization period. The state of the building zone includes atemperature of air within the building zone and a temperature of solidmass within the building zone. The method includes using a predictivemodel of the building zone to predict the state of the building zone ateach of a plurality of time steps of the optimization period based onthe estimated initial state of the building zone and a temperaturesetpoint trajectory including a temperature setpoint for each of theplurality of time steps. The method includes generating optimaltemperature setpoints for each of the plurality of time steps byoptimizing a cost function that accounts for a cost operating HVACequipment during each of the plurality of time steps.

In some embodiments, the method includes generating the cost functionand determining the cost operating the HVAC equipment during each timestep using a set of time-varying utility rates comprising a utility ratevalue for each time step. The set of time-varying utility rates can bereceived from a utility provider or predicted as part of the method.

In some embodiments, the predicted state of the building zone includes apredicted temperature of the building zone. The method may includeoptimizing the cost function subject to a constraint on the predictedtemperature of the building zone.

In some embodiments, the method includes predicting a value of anunmeasured heat load disturbance experienced by the building zone ateach of the plurality of time steps in the optimization period. In someembodiments, the predictive model defines the temperature of thebuilding zone as a function of the unmeasured heat load disturbanceexperienced by the building zone at each time step.

In some embodiments, the predictive model includes a thermal massstorage model that defines the temperature of the building zone as afunction of heat transfer between air within the building zone and solidmass within the building zone. In some embodiments, the predictive modeldefines the temperature of the building zone as a function of heattransfer between the building zone and the HVAC equipment. In someembodiments, the predictive model includes an HVAC load model thatdefines an amount of heating or cooling provided by the HVAC equipmentas a function of the temperature of the building zone and thetemperature setpoint.

In some embodiments, the method includes predicting the cost ofoperating the HVAC equipment as a function of an amount of heating orcooling provided by the HVAC equipment at each time step of theoptimization period.

In some embodiments, the method includes generating the predictive modelby performing a system identification process. The system identificationprocess includes collecting a set of input-output data and fittingparameters of the predictive model to the set of input-output data. Theinput-output data may include values of the temperature setpoint and thetemperature of the building zone during each of a plurality of timesteps during a learning period.

Another implementation of the present disclosure is a model predictivecontroller for monitoring and controlling temperature of a buildingzone. The model predictive controller includes a state/disturbanceestimator and a predictive optimizer. The state/disturbance estimator isconfigured to estimate an initial state of the building zone at abeginning of an optimization period. The state of the building zoneincludes a temperature of air within the building zone and a temperatureof solid mass within the building zone. The predictive optimizer isconfigured to use a predictive model of the building zone to predict thestate of the building zone at each of a plurality of time steps of theoptimization period based on the estimated initial state of the buildingzone and a temperature setpoint trajectory including a temperaturesetpoint for each of the plurality of time steps. The predictiveoptimizer is configured to generate temperature setpoints for each ofthe plurality of time steps by optimizing a cost function that accountsfor a cost operating HVAC equipment during each of the plurality of timesteps.

Another implementation of the present disclosure is a thermostat formonitoring and controlling temperature of a building zone. Thethermostat includes a load/rate predictor and a predictive optimizer.The load/rate predictor is configured to predict a price of one or moreresources consumed by HVAC equipment to generate heating or cooling forthe building zone at each of a plurality of time steps of anoptimization period. The predictive optimizer is configured to generatea cost function that accounts for a cost of operating the HVAC equipmentduring each time step as a function of the predicted prices at each timestep, use a predictive model to predict the temperature of the buildingzone during each of the plurality of time steps as a function of atemperature setpoint trajectory including a temperature setpoint foreach of the plurality of time steps, and optimize the cost functionsubject to a constraint on the predicted temperature of the buildingzone to determine optimal temperature setpoints for each of theplurality of time steps.

In some embodiments, the thermostat includes an equipment controllerconfigured to drive the temperature of the building zone to an optimaltemperature setpoint by operating the HVAC equipment to provide heatingor cooling to the building zone.

In some embodiments, the predictive model defines the temperature of thebuilding zone as a function of an unmeasured heat load disturbance. Insome embodiments, the load/rate predictor is configured to predict avalue of the unmeasured heat load disturbance experienced by thebuilding zone at each of the plurality of time steps in the optimizationperiod as a function of at least one of day type, time of day, buildingoccupancy, outdoor air temperature, and weather forecasts.

In some embodiments, the thermostat includes a system identifierconfigured to generate the predictive model by performing a systemidentification process. The system identification process includescollecting a set of input-output data and fitting parameters of thepredictive model to the set of input-output data. The input-output datamay include values of the temperature setpoint and the temperature ofthe building zone during each of a plurality of time steps during alearning period.

In some embodiments, the system identification process includesmodulating the temperature setpoint within a constrained temperaturesetpoint range and recording values of the modulated temperaturesetpoint and values of the temperature of the building zone that resultfrom modulating the temperature setpoint.

In some embodiments, the predictive model includes a thermal massstorage model that defines the temperature of the building zone as afunction of heat transfer between air within the building zone and solidmass within the building zone. In some embodiments, the predictive modeldefines the temperature of the building zone as a function of heattransfer between the building zone and the HVAC equipment. In someembodiments, the predictive model includes an HVAC load model thatdefines the heating or cooling provided by the HVAC equipment as afunction of the temperature of the building zone and the temperaturesetpoint.

In some embodiments, the predictive optimizer is configured to predictthe cost of operating the HVAC equipment as a function of an amount ofthe heating or cooling provided by the HVAC equipment at each time stepof the optimization period.

Another implementation of the present disclosure is a method performedby a thermostat for a building zone for monitoring and controllingtemperature of the building zone. The method includes predicting a priceof one or more resources consumed by HVAC equipment to generate heatingor cooling for the building zone at each of a plurality of time steps ofan optimization period, generating a cost function that accounts for acost of operating the HVAC equipment during each time step as a functionof the predicted prices at each time step, using a predictive model topredict the temperature of the building zone during each of theplurality of time steps as a function of a temperature setpointtrajectory including a temperature setpoint for each of the plurality oftime steps, and optimizing the cost function subject to a constraint onthe predicted temperature of the building zone to determine optimaltemperature setpoints for each of the plurality of time steps.

In some embodiments, the method includes driving the temperature of thebuilding zone to an optimal temperature setpoint by operating the HVACequipment to provide heating or cooling to the building zone.

In some embodiments, the predictive model defines the temperature of thebuilding zone as a function of an unmeasured heat load disturbance. Insome embodiments, the method includes predicting a value of theunmeasured heat load disturbance experienced by the building zone ateach of the plurality of time steps in the optimization period as afunction of at least one of day type, time of day, building occupancy,outdoor air temperature, and weather forecasts.

In some embodiments, the method includes generating the predictive modelby performing a system identification process. The system identificationprocess may include collecting a set of input-output data and fittingparameters of the predictive model to the set of input-output data. Theinput-output data may include values of the temperature setpoint and thetemperature of the building zone during each of a plurality of timesteps during a learning period.

In some embodiments, performing the system identification processincludes modulating the temperature setpoint within a constrainedtemperature setpoint range and recording values of the modulatedtemperature setpoint and values of the temperature of the building zonethat result from modulating the temperature setpoint.

In some embodiments, the predictive model includes a thermal massstorage model that defines the temperature of the building zone as afunction of heat transfer between air within the building zone and solidmass within the building zone. In some embodiments, the predictive modeldefines the temperature of the building zone as a function of heattransfer between the building zone and the HVAC equipment. In someembodiments, the predictive model includes an HVAC load model thatdefines the heating or cooling provided by the HVAC equipment as afunction of the temperature of the building zone and the temperaturesetpoint.

Another implementation of the present disclosure is a thermostat formonitoring and controlling temperature of a building zone. Thethermostat includes a load/rate predictor and a predictive optimizer.The load/rate predictor is configured to predict a price of one or moreresources consumed by HVAC equipment to generate heating or cooling forthe building zone at each of a plurality of time steps of anoptimization period. The predictive optimizer is configured to generatea cost function that accounts for a cost of operating the HVAC equipmentduring each time step as a function of the predicted prices at each timestep, use a predictive model to predict the temperature of the buildingzone during each of the plurality of time steps as a function of atemperature setpoint trajectory including a temperature setpoint foreach of the plurality of time steps, and perform an optimization of thecost function subject to a constraint on the predicted temperature ofthe building zone to determine temperature setpoints for each of theplurality of time steps.

Those skilled in the art will appreciate that the summary isillustrative only and is not intended to be in any way limiting. Otheraspects, inventive features, and advantages of the devices and/orprocesses described herein, as defined solely by the claims, will becomeapparent in the detailed description set forth herein and taken inconjunction with the accompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a drawing of a smart thermostat having a landscape aspectratio, according to some embodiments.

FIG. 2 is a drawing of a smart thermostat having a portrait aspectratio, according to some embodiments.

FIG. 3 is a block diagram of a model predictive control system includinga model predictive controller and a smart thermostat which communicateswith the model predictive controller via a network, according to someembodiments.

FIG. 4 is a block diagram of a model predictive control system includinga smart thermostat that includes a model predictive controller and anequipment controller, according to some embodiments.

FIG. 5 is a block diagram illustrating the model predictive controllerof FIGS. 3-4 in greater detail, according to some embodiments.

FIG. 6 is a block diagram illustrating various types of predictivemodels, according to some embodiments.

FIG. 7 is a graph of a temperature excitation signal over a one-dayperiod, according to some embodiments.

FIG. 8 is a thermal circuit diagram that models the heat transfercharacteristics of a building zone, according to some embodiments.

FIG. 9 is a flowchart of a process which can be performed by the modelpredictive controller of FIG. 5 to identify parameters of a predictivemodel, according to some embodiments.

FIG. 10 is a graph of an unfiltered temperature signal and an unfilteredsensible load signal which can be received as feedback from a controlledHVAC system, according to some embodiments.

FIG. 11 is a graph of a filtered temperature signal and a filteredsensible load signal which can be generated by filtering the signalsshown in FIG. 10, according to some embodiments.

FIG. 12 is a flowchart of a process for operating the model predictivecontroller of FIG. 5 in a system identification mode and in anoperational mode, according to some embodiments.

FIG. 13 is a flowchart of a process which can be performed by the modelpredictive controller of FIG. 5 to determine optimal temperaturesetpoints when operating in the operational mode, according to someembodiments.

FIG. 14 is a flowchart of a process which can be performed by the modelpredictive controller of FIG. 5 to generate a predictive model and usethe predictive model to determine optimal temperature setpoints for abuilding zone, according to some embodiments.

FIG. 15 is a graph comparing the energy consumption of the modelpredictive control systems of FIGS. 3-4 with a baseline system that doesnot use model predictive control, according to some embodiments.

FIG. 16 is a graph of time-varying electricity prices illustratingdifferent electricity prices at different times, according to someembodiments.

FIG. 17 is a graph comparing the building mass temperature in the modelpredictive control systems of FIGS. 3-4 with a baseline system that doesnot use model predictive control, according to some embodiments.

FIG. 18 is a graph comparing the zone air temperature in the modelpredictive control systems of FIGS. 3-4 with a baseline system that doesnot use model predictive control, according to some embodiments.

DETAILED DESCRIPTION Overview

Referring generally to the FIGURES, a smart thermostat with modelpredictive control and components thereof are shown according, to someembodiments. The smart thermostat can be configured to monitor andcontrol one or more environmental conditions of a building zone (e.g.,temperature, humidity, air quality, etc.). In some embodiments, thesmart thermostat is mounted on a wall within the building zone andconfigured to measure the temperature, humidity, and/or otherenvironmental conditions of the building zone. The smart thermostat canbe configured to communicate with HVAC equipment that operates to affectthe measured environmental conditions.

In some embodiments, the smart thermostat includes a model predictivecontroller and an equipment controller. In other embodiments, the modelpredictive controller is separate from the smart thermostat andcommunicates with the smart thermostat via a communications network(e.g., the Internet, a building network, etc.). The model predictivecontroller can be configured to determine optimal temperature setpointsT_(sp) for the building zone for each of a plurality of time steps in anoptimization period. The equipment controller can be configured toreceive the optimal temperature setpoints T_(sp) from the modelpredictive controller and can operate HVAC equipment to drive thetemperature of the building zone to the optimal temperature setpoints.

To determine the optimal temperature setpoints T_(sp), the modelpredictive controller can optimize an objective function (i.e., a costfunction) that accounts for the cost of operating the HVAC equipmentover the duration of the optimization period. The costs of operating theHVAC equipment can include, for example, the costs of resources consumedby the HVAC equipment during operation (e.g., electricity, natural gas,water, etc.), demand charges imposed by an electric utility, peak loadcontribution charges, equipment degradation/replacement costs, and/orother costs associated with the operation of the HVAC equipment. Theoptimization performed by the model predictive controller is describedin greater detail below.

The model predictive controller can optimize the objective functionsubject to a set of constraints. The constraints may include temperatureconstraints (e.g., a maximum temperature limit and a minimum temperaturelimit for the building zone), equipment capacity constraints, loadchange constraints, thermal mass storage constraints, HVAC loadconstraints, and/or other constraints that limit the operation of theHVAC equipment and/or describe the temperature evolution of the buildingzone. The model predictive controller can automatically generate theoptimization constraints by performing a system identification processand generating predictive models that describe the controlled system(i.e., the HVAC equipment and the building zone).

In some embodiments, the model predictive controller generates a thermalmass storage model that describes the temperature of the air within thebuilding zone T_(ia) as a function of several physical parameters andvariables. The physical parameters in the thermal mass storage model mayinclude, for example, the thermal capacitance C_(ia) of the air in thebuilding zone, the thermal capacitance C_(m) of the solid mass withinthe building zone, the thermal resistance R_(mi) between the air withinthe building zone and the solid mass, the thermal resistance R_(oi)between the air within the building zone and the ambient environment,etc. The variables in the thermal mass storage model may include thetemperature of air within the building zone T_(ia), the temperature ofthe solid mass within the building zone T_(m), the outside airtemperature T_(oa), the amount of heating or cooling provided by theHVAC equipment {dot over (Q)}_(HVAC), and/or a thermal energy loaddisturbance (i.e., a heat load disturbance) {dot over (Q)}_(other).

In some embodiments, the model predictive controller generates an HVACload model. The HVAC load model may describe the amount of heating orcooling {dot over (Q)}_(HVAC) provided by the HVAC equipment as afunction of the temperature setpoint T_(sp) and the temperature of theair T_(ia) within the building zone. The model predictive controller canuse the thermal mass storage model and the HVAC load model to predictsystem states (i.e., building mass temperature T_(m), building airtemperature T_(ia), etc.), predict the load disturbance {dot over(Q)}_(other), and establish constraints on the optimization performed bythe model predictive controller.

In some embodiments, the model predictive controller uses the thermalmass storage model and the HVAC load model to take advantage of thebuilding mass as a thermal energy storage medium. For example, thermalenergy can be stored in the solid building mass (i.e., precooling orpreheating the building zone) to allow the model predictive controllerto take advantage of time-varying utility rates and demand charges whendetermining the optimal temperature setpoints T_(sp). To precool thebuilding zone, thermal energy may be removed from the thermal massduring periods when energy prices are lowest. During periods when energyprices are largest, thermal energy may be moved back into the solid massfrom the zone air, which reduces the sensible HVAC load {dot over(Q)}_(HVAC) needed to maintain the zone air within a comfortabletemperature range. Similarly, the model predictive controller canpreheat the solid mass within the building zone when energy prices arelow and use the stored thermal energy to warm the zone air when energyprices are higher.

In some embodiments, the model predictive controller is configured tooperate in two distinct modes: (1) parameter identification mode and (2)operational mode. In the parameter identification mode, the modelpredictive controller can manipulate the temperature setpoint T_(sp)provided to the smart thermostat and/or the equipment controller toinduce a dynamic response in the zone temperature T_(ia). This procedureof inducing a dynamic response is referred to as the parameteridentification experiment or system identification process and can beused to identify the values of the parameters in the thermal massstorage model and the HVAC load model. In the operational mode, themodel predictive controller can solve an optimal control problem todetermine the optimal temperature setpoint trajectory T_(sp) (i.e., atime series of temperature setpoints) that minimizes the cost of energyconsumed by the HVAC equipment over the duration of the optimizationperiod. The optimal temperature setpoints T_(sp) can be provided to thesmart thermostat and/or the equipment controller for use in generatingoperating commands for the HVAC equipment.

The setpoint T_(sp) computed for the first time step of the optimizationperiod can be sent to the smart thermostat and/or the equipmentcontroller to be implemented during the first time step. The smartthermostat and/or the equipment controller can turn stages of HVACequipment on and/or off such that the temperature T_(ia) of the zone isforced to and then maintained at the temperature setpoint T_(sp). At thenext time step, the model predictive controller receives updatedinformation (i.e., feedback on the zone conditions), resolves theoptimal control problem, and sends the temperature setpoint T_(sp) forthe next time step to the smart thermostat and/or the equipmentcontroller.

The model predictive controller can collect feedback on the zoneconditions (e.g., zone temperature T_(ia), outdoor air temperatureT_(oa), heating/cooling load {dot over (Q)}_(HVAC)) when operating inthe parameter identification mode. The zone conditions can be monitoredand recorded over the identification experiment, along with thecommanded zone temperature setpoint T_(sp), as input-output data (i.e.,training data). The model predictive controller can use the input-outputdata to compute the zone thermal model parameters by performing a systemidentification process. The feedback received from the smart thermostatmay include the current minimum and maximum allowable setpoint to limitthe adjustments to the zone temperature setpoint T_(sp) made during theidentification experiment. These and other features of the smartthermostat and/or the model predictive controller are described ingreater detail below.

Throughout this disclosure, the terms “optimal,” “optimized,” “optimum,”and the like are used to refer to values (e.g., temperature setpointvalues, HVAC load values, etc.) that are determined by performing anoptimization process. Similarly, the term “optimizing” is used to referto the process of performing an optimization. In some instances, thevalues resulting from the optimization process are true optimal values(i.e., values that achieve the minimum possible value or maximumpossible value for a performance variable given the constraints on theoptimization process). In other instances, the values resulting from theoptimization process are not true optimal values. This can result fromimperfect information used to perform the optimization, inaccuracies inthe predictive model used to constrain the optimization, and/or variousother factors that can prevent the optimization from converging on thetrue optimal values. The terms “optimal,” “optimized,” “optimum,” andthe like should be interpreted to include any values that result fromperforming an optimization, regardless of whether the values are trueoptimal values. Similarly, the term “optimizing” should be interpretedto include the process of performing an optimization, regardless ofwhether that optimization converges on the true optimal values.

Smart Thermostat

Referring now to FIGS. 1-2, a smart thermostat 100 is shown, accordingto some embodiments. Smart thermostat 100 can be configured to monitorand control one or more environmental conditions of a building zone(e.g., temperature, humidity, air quality, etc.). In some embodiments,smart thermostat 100 is mounted on a wall within the building zone andconfigured to measure the temperature, humidity, and/or otherenvironmental conditions of the building zone. Smart thermostat 100 canbe configured to communicate with HVAC equipment that operates to affectthe measured environmental conditions.

Smart thermostat 100 can be installed in a home, office building,school, hospital, or any other environment-controlled space. Inresidential implementations, the HVAC equipment controlled by smartthermostat 100 may include a home furnace, air conditioning unit, and/orother types of residential HVAC equipment. In commercialimplementations, the HVAC equipment may include one or more chillers,boilers, air handling units, rooftop units, dampers, or other types ofHVAC equipment configured to affect the environment of a building zone.It is contemplated that smart thermostat 100 can be configured tocontrol any type of HVAC equipment.

In some embodiments, smart thermostat 100 includes a communicationsinterface that enables smart thermostat 100 to connect to acommunications network (e.g., a local area network, the Internet, acellular network, etc.) and communicate with various external systemsand devices (e.g., user devices, remote servers, remote controllers,etc.). For example, smart thermostat 100 can be configured to receiveweather forecasts from a remote weather service via the Internet. Insome embodiments, smart thermostat 100 receives temperature setpointsfrom a model predictive controller via the communications network. Inother embodiments, the model predictive controller is a component ofsmart thermostat 100.

Smart thermostat 100 is shown to include a display screen 102 and a base104. Base 104 can attach to a wall or other surface upon which smartthermostat 100 is mounted. Display screen 102 may be transparent orsemi-transparent (e.g., an organic LED display) and can be configured todisplay text, graphics, and other information for presentation to auser. In some embodiments, display screen 102 is touch-sensitive (e.g.,a capacitive or resistive touch screen) and configured to receive userinput. Display screen 102 may have a landscape aspect ratio as shown inFIG. 1 (i.e., the width of display screen 102 may exceed the height ofdisplay screen 102) or a portrait aspect ratio as shown in FIG. 2 (i.e.,the height of display screen 102 may exceed the width of display screen102). Display screen 102 may be cantilevered from base 104 (i.e.,attached to base 104 along a single edge of display screen 102) and mayextend from base 104 in a direction substantially parallel to themounting surface.

Smart thermostat 100 can include some of all of the features of thethermostats described in U.S. patent application Ser. No. 15/143,373filed Apr. 29, 2016, U.S. patent application Ser. No. 15/146,763 filedMay 4, 2016, U.S. patent application Ser. No. 15/146,749 filed May 4,2016, U.S. patent application Ser. No. 15/146,202 filed May 4, 2016,U.S. patent application Ser. No. 15/146,134 filed May 4, 2016, U.S.patent application Ser. No. 15/146,649 filed May 4, 2016, U.S.Provisional Patent Application No. 62/331,863 filed May 4, 2016, U.S.Provisional Patent Application No. 62/352,955 filed Jun. 21, 2016, U.S.patent application Ser. No. 15/298,191 filed Oct. 19, 2016, U.S. patentapplication Ser. No. 15/336,793 filed Oct. 28, 2016, U.S. patentapplication Ser. No. 15/336,792 filed Oct. 28, 2016, U.S. patentapplication Ser. No. 15/336,789 filed Oct. 28, 2016, U.S. patentapplication Ser. No. 15/338,221 filed Oct. 28, 2016, U.S. patentapplication Ser. No. 15/338,215 filed Oct. 28, 2016, U.S. patentapplication Ser. No. 15/336,791 filed Oct. 28, 2016, U.S. patentapplication Ser. No. 15/397,722 filed Jan. 3, 2017, and/or U.S.Provisional Patent Application No. 62/446,296 filed Jan. 13, 2017. Theentire disclosure of each of these patent applications is incorporatedby reference herein.

Model Predictive Control Systems

Referring now to FIGS. 3-4, a pair of model predictive control (MPC)systems 300 and 400 are shown, according to some embodiments. MPC system300 is shown to include a model predictive controller 302, acommunications network 304, smart thermostat 100, HVAC equipment 308,and a building zone 310. In MPC system 300, model predictive controller302 is separate from smart thermostat 100 and configured to provide atemperature setpoint T_(sp) to smart thermostat 100 via communicationsnetwork 304. In MPC system 400, model predictive controller 302 is acomponent of smart thermostat 100 along with an equipment controller406.

Building zone 310 can include one or more rooms or zones within a home,office building, school, hospital, or any other environment-controlledspace. HVAC equipment 308 can include any type of equipment operable toaffect the temperature, humidity, and/or other environmental conditionsof building zone 310. For example, HVAC equipment 308 can include a homefurnace, air conditioning unit, one or more chillers, boilers, airhandling units, rooftop units, dampers, or other types of HVAC equipmentconfigured to affect the environment of building zone 310.

Model predictive controller 302 can be configured to determine anoptimal temperature setpoint T_(sp) for smart thermostat 100 for each ofa plurality of time steps during an optimization period. Smartthermostat 100 can use the temperature setpoints T_(sp) provided bymodel predictive controller 302 to generate equipment commands for HVACequipment 308. The equipment commands can be generated by smartthermostat 100 and/or equipment controller 406 (e.g., an on/offcontroller, an equipment staging controller, a proportional-integral(PI) controller, a proportional-integral-derivative (PID) controller,etc.). HVAC equipment 308 operate according to the equipment commands toprovide variable amount of heating or cooling {dot over (Q)}_(HVAC) tobuilding zone 310. By controlling the temperature setpoint T_(sp), modelpredictive controller 302 can modulate the amount of heating or cooling{dot over (Q)}_(HVAC) provided by HVAC equipment 308, thereby affectingthe temperature of building zone 310.

To determine the optimal temperature setpoints T_(sp), model predictivecontroller 302 can optimize an objective function (i.e., a costfunction) that accounts for the cost of operating HVAC equipment 308over the duration of the optimization period. The costs of operatingHVAC equipment 308 can include, for example, the costs of resourcesconsumed by HVAC equipment 308 during operation (e.g., electricity,natural gas, water, etc.), demand charges imposed by an electricutility, peak load contribution charges, equipmentdegradation/replacement costs, and/or other costs associated with theoperation of HVAC equipment 308. The optimization performed by modelpredictive controller 302 is described in greater detail with referenceto FIG. 5.

Model predictive controller 302 can optimize the objective functionsubject to a set of constraints. The constraints may include temperatureconstraints (e.g., a maximum temperature limit and a minimum temperaturelimit for building zone 310), equipment capacity constraints, loadchange constraints, thermal mass storage constraints, HVAC loadconstraints, and/or other constraints that limit the operation of HVACequipment 308 and/or describe the temperature evolution of building zone310. Model predictive controller 302 can automatically generate theoptimization constraints by performing a system identification processand generating predictive models that describe the controlled system(i.e., HVAC equipment 308 and building zone 310).

In some embodiments, model predictive controller 302 generates a thermalmass storage model that describes the temperature of the air withinbuilding zone 310 T_(ia) as a function of several physical parametersand variables. The physical parameters in the thermal mass storage modelmay include, for example, the thermal capacitance C_(ia) of the air inbuilding zone 310, the thermal capacitance C_(m) of the solid masswithin building zone 310, the thermal resistance R_(mi) between the airwithin building zone 310 and the solid mass, the thermal resistanceR_(oi) between the air within building zone 310 and the ambientenvironment, etc. The variables in the thermal mass storage model mayinclude the temperature of air within building zone 310 T_(ia), thetemperature of the solid mass within building zone 310 T_(m), theoutside air temperature T_(oa), the amount of heating or coolingprovided by HVAC equipment 308 {dot over (Q)}_(HVAC), and/or a thermalenergy load disturbance {dot over (Q)}_(other).

In some embodiments, model predictive controller 302 generates a HVACload model. The HVAC load model may describe the amount of heating orcooling {dot over (Q)}_(HVAC) provided by HVAC equipment 308 as afunction of the temperature setpoint T_(sp) and the temperature of theair T_(ia) within building zone 310. Model predictive controller 302 canuse the thermal mass storage model and the HVAC load model to predictsystem states (i.e., building mass temperature T_(m), building airtemperature T_(ia), etc.), predict the load disturbance {dot over(Q)}_(other), and establish constraints on the optimization performed bymodel predictive controller 302.

In some embodiments, model predictive controller 302 uses the thermalmass storage model and the HVAC load model to take advantage of thebuilding mass as a thermal energy storage medium. For example, thermalenergy can be stored in the solid building mass (i.e., precooling orpreheating building zone 310) to allow model predictive controller 302to take advantage of time-varying utility rates and demand charges whendetermining the optimal temperature setpoints T_(sp). To precoolbuilding zone 310, thermal energy may be removed from the thermal massduring periods when energy prices are lowest. During periods when energyprices are largest, thermal energy may be moved back into the solid massfrom the zone air, which reduces the sensible HVAC load {dot over(Q)}_(HVAC) needed to maintain the zone air within a comfortabletemperature range. Similarly, model predictive controller 302 canpreheat the solid mass within building zone 310 when energy prices arelow and use the stored thermal energy to warm the zone air when energyprices are higher.

In some embodiments, model predictive controller 302 is configured tooperate in two distinct modes: (1) parameter identification mode and (2)operational mode. In the parameter identification mode, model predictivecontroller 302 can manipulate the temperature setpoint T_(sp) providedto smart thermostat 100 and/or equipment controller 406 to induce adynamic response in the zone temperature T_(ia). This procedure ofinducing a dynamic response is referred to as the parameteridentification experiment or system identification process and can beused to identify the values of the parameters in the thermal massstorage model and the HVAC load model. In the operational mode, modelpredictive controller 302 can solve an optimal control problem todetermine the optimal temperature setpoint trajectory T_(sp) (i.e., atime series of temperature setpoints) that minimizes the cost of energyconsumed by HVAC equipment 308 over the duration of the optimizationperiod. The optimal temperature setpoints T_(sp) can be provided tosmart thermostat 100 and/or equipment controller 406 for use ingenerating operating commands for HVAC equipment 308.

The setpoint T_(sp) computed for the first time step of the optimizationperiod can be sent to smart thermostat 100 and/or equipment controller406 to be implemented during the first time step. Smart thermostat 100and/or equipment controller 406 can turn stages of HVAC equipment 308 onand/or off such that the temperature T_(ia) of the zone is forced to andthen maintained at the temperature setpoint T_(sp). At the next timestep, model predictive controller 302 receives updated information(i.e., feedback on the zone conditions), resolves the optimal controlproblem, and sends the temperature setpoint T_(sp) for the next timestep to smart thermostat 100 and/or equipment controller 406.

Model predictive controller 302 can collect feedback on the zoneconditions (e.g., zone temperature T_(ia), outdoor air temperatureT_(oa), heating/cooling load {dot over (Q)}_(HVAC)) when operating inthe parameter identification mode. The zone conditions can be monitoredand recorded over the identification experiment, along with thecommanded zone temperature setpoint T_(sp), as input-output data (i.e.,training data). Model predictive controller 302 can use the input-outputdata to compute the zone thermal model parameters by performing a systemidentification process. The feedback received from smart thermostat 100may include the current minimum and maximum allowable setpoint to limitthe adjustments to the zone temperature setpoint T_(sp) made during theidentification experiment.

In some embodiments, a pseudo-random binary signal (PRBS) is generatedthat is subsequently used to produce a signal that takes values in theset {T_(sp,1), T_(sp,2)}, where T_(sp,1) and T_(sp,2) are maximum andminimum allowable temperature setpoints depending on mode of operationof smart thermostat 100 (e.g., heating, cooling, home, sleep, or awaymode). The minimum and maximum temperature setpoints may be time-varyingto account for the different operational modes of smart thermostat 100(e.g., home, sleep, and away). In some embodiments, smart thermostat 100uses occupancy detection and/or user feedback to change the mode.

To determine the optimal temperature setpoints T_(sp) in the operationalmode, model predictive controller 302 can use the thermal model ofbuilding zone 310 to predict the zone temperature T_(ia) and HVAC fuelconsumption over the optimization period, given a forecast of theweather and heat disturbance load on the zone (e.g., heat generated frompeople and electrical equipment and gained through solar radiation). Thezone temperature T_(ia) may be subject to comfort constraints. Thepredicted trajectories of the zone temperature T_(ia), temperaturesetpoint T_(sp), and HVAC fuel consumption can be sent to smartthermostat 100 to be displayed to the user allowing the user to betterunderstand the setpoint decision made by model predictive controller302.

In some embodiments, smart thermostat 100 allows a user to override theoptimal temperature setpoints computed by model predictive controller302. Smart thermostat 100 may include a savings estimator configured todetermine a potential cost savings resulting from the optimaltemperature setpoints relative to the user-specified temperaturesetpoints. Smart thermostat 100 can be configured to display thepotential cost savings via display screen 102 to inform the user of theeconomic cost predicted to result from overriding the optimaltemperature setpoints. These and other features of model predictivecontroller 302 are described in greater detail below.

Model Predictive Controller

Referring now to FIG. 5, a block diagram illustrating model predictivecontroller 302 in greater detail is shown, according to comeembodiments. Model predictive controller 302 is shown to include acommunications interface 502 and a processing circuit 504.Communications interface 502 may facilitate communications between modelpredictive controller 302 and external systems or devices. For example,communications interface 502 may receive measurements of zone conditionsfrom smart thermostat 100, HVAC equipment 308, and/or building zone 310.Zone conditions may include the temperature T_(ia) of the air withinbuilding zone 310, the outside air temperature T_(oa), and/or the amountof heating/cooling {dot over (Q)}_(HVAC) provided to building zone 310at each time step of the optimization period. Communications interface502 may also receive temperature constraints for building zone 310and/or equipment constraints for HVAC equipment 308. In someembodiments, communications interface 502 receives utility rates fromutilities 524 (e.g., time-varying energy prices, demand charges, etc.),weather forecasts from a weather service 526, and/or historical load andrate data from historical data 528. In some embodiments, modelpredictive controller 302 uses communications interface 502 to providetemperature setpoints T_(sp) to smart thermostat 100 and/or equipmentcontroller 406.

Communications interface 502 may include wired or wirelesscommunications interfaces (e.g., jacks, antennas, transmitters,receivers, transceivers, wire terminals, etc.) for conducting datacommunications external systems or devices. In various embodiments, thecommunications may be direct (e.g., local wired or wirelesscommunications) or via a communications network (e.g., a WAN, theInternet, a cellular network, etc.). For example, communicationsinterface 502 can include an Ethernet card and port for sending andreceiving data via an Ethernet-based communications link or network. Inanother example, communications interface 502 can include a Wi-Fitransceiver for communicating via a wireless communications network orcellular or mobile phone communications transceivers.

Processing circuit 504 is shown to include a processor 506 and memory508. Processor 506 may be a general purpose or specific purposeprocessor, an application specific integrated circuit (ASIC), one ormore field programmable gate arrays (FPGAs), a group of processingcomponents, or other suitable processing components. Processor 506 isconfigured to execute computer code or instructions stored in memory 508or received from other computer readable media (e.g., CDROM, networkstorage, a remote server, etc.).

Memory 508 may include one or more devices (e.g., memory units, memorydevices, storage devices, etc.) for storing data and/or computer codefor completing and/or facilitating the various processes described inthe present disclosure. Memory 508 may include random access memory(RAM), read-only memory (ROM), hard drive storage, temporary storage,non-volatile memory, flash memory, optical memory, or any other suitablememory for storing software objects and/or computer instructions. Memory508 may include database components, object code components, scriptcomponents, or any other type of information structure for supportingthe various activities and information structures described in thepresent disclosure. Memory 508 may be communicably connected toprocessor 506 via processing circuit 504 and may include computer codefor executing (e.g., by processor 506) one or more processes describedherein. When processor 506 executes instructions stored in memory 508for completing the various activities described herein, processor 506generally configures controller 302 (and more particularly processingcircuit 504) to complete such activities.

Still referring to FIG. 5, model predictive controller 302 is shown toinclude a system identifier 510, a load/rate predictor 518, astate/disturbance estimator 520, and a predictive optimizer 522. Inbrief overview, system identifier 510 can be configured to perform asystem identification process to identify the values of the parametersin the building thermal mass model and the HVAC load model. Theparameters identified by system identifier 510 may define the values ofthe system matrices A, B, C, and D and the estimator gain K.State/disturbance estimator 520 can be configured to compute an estimateof the current state {circumflex over (x)}(k) and unmeasured disturbance{circumflex over (d)}(k). The estimation performed by state/disturbanceestimator 520 can be performed at each sample time and may use thesystem matrices A, B, C, D and estimator gain K identified by systemidentifier 510.

Load/rate predictor 518 can generate predictions or forecasts of theinternal heat load {dot over (Q)}_(other), the non-HVAC buildingelectricity consumption (used to calculate electrical demand charges),and the utility rates at each time step of the optimization period as afunction of historical data 528 and several key predictor variables(e.g., time-of-day, outdoor air conditions, etc.). The models used byload/rate predictor 518 may be different than that identified by systemidentifier 510. In some embodiments, load/rate predictor 518 uses theprediction techniques described in U.S. patent application Ser. No.14/717,593 filed May 20, 2015, the entire disclosure of which isincorporated by reference herein.

Predictive optimizer 522 can manipulate the temperature setpoint T_(sp)of building zone 310 to minimize the economic cost of operating HVACequipment 308 over the duration of the optimization period using thebuilding thermal mass storage model, the HVAC load model, and state anddisturbance estimates and forecasts. In some embodiments, predictiveoptimizer 522 uses open-loop predictions to optimize the temperaturesetpoint T_(sp) trajectory over the optimization period. At each timestep of the optimization period, the predictions can be updated based onfeedback from the controlled system (e.g., measured temperatures,measured HVAC loads, etc.).

Before discussing system identifier 510, load/rate predictor 518,state/disturbance estimator 520, and predictive optimizer 522 in detail,the notation used throughout the remainder of the present disclosure andthe class of system model used by model predictive controller 302 areexplained. The set of integers is denoted by

and the set of positive integers is denoted by

_(≥0), and the set of integers contained in the interval [a, b] isdenoted by

_(a:b). For a time-dependent vector x(k)∈

^(n), {circumflex over (x)}(i|k)∈

^(n) denotes the estimated value of x at time step i∈

_(≥0) given the measurement at time step k∈

_(≥0) where i≥k and {tilde over (x)}(i|k)∈

^(n) denotes the (open-loop) predicted value of x at time step i∈

_(≥0) with the prediction initialized at time k∈

_(≥0) (i≥k). For notational simplicity, the notation of the predictedvalue of x at time i starting from time k is abbreviated to {tilde over(x)}(i). Boldface letters are used to represent a sequence withcardinality N∈

_(≥0) (i.e., x:={x(0), . . . , x(N−1)}). The notation ⋅* (e.g., x*)denotes an optimal quantity with respect to some optimization problem.

In some embodiments, the model generated and used by model predictivecontroller 302 is a discrete-time linear time-invariant model, as shownin the following equation:

x(k+1)=Ax(k)+Bu(k)

y(k)=Cx(k)+Du(k)  (Equation 1)

where k∈

_(≥0) is the time index, x(k)∈

^(n) is the state vector, y(k)∈

^(p) is the measured output, and u(k)∈

^(m) is the input vector. The matrix A can be assumed to be stable; thatis the real parts of the eigenvalues of A lie within the unit circle.The pair (A, B) can be assumed to be controllable and the pair (A, C)can be assumed to be observable. Some HVAC equipment 308 and devices areonly operated at finite operating stages. In the building modelsconsidered, this gives rise to discrete outputs whereby the measuredoutput takes values within a finite set of numbers:y_(i)(k)∈{y_(i,stage,1), y_(i,stage,2), . . . , y_(i,stage,m)} for somei. In this case, the output does not evolve according to Equation 1.

In some embodiments, the model generated and used by model predictivecontroller 302 is a grey box continuous-time models derived from thermalresistance-capacitance (RC) modeling principals. Model predictivecontroller 302 can be configured to parameterize the resultingcontinuous-time thermal models and identify the continuous-time modelparameters. The continuous-time version of the model of Equation 1 canbe represented as shown in the following equation

{dot over (x)}(t)=A _(c) x(t)+B _(c) u(t)

y(t)=C _(c) x(t)+D _(c) u(t)  (Equation 2)

where t≥0 is continuous time (the initial time is taken to be zero) andA_(c), B_(c), C_(c), D_(c) are the continuous-time versions of thediscrete-time system matrices A, B, C, and D. The values of theparameters in the matrices A, B, C, and D as well as an estimator gain Kcan be determined by system identifier 510 by performing a systemidentification process (described in greater detail below)

State/Disturbance Estimator

State/disturbance estimator 520 can be configured to compute an estimateof the current state {circumflex over (x)}(k) and unmeasured disturbance{circumflex over (d)}(k). The estimation performed by state/disturbanceestimator 520 can be performed at each sample time (i.e., when operatingin the operational mode) and may use the system matrices A, B, C, D andestimator gain K identified by system identifier 510. In someembodiments, state/disturbance estimator 520 modifies the system modelto account for measurement noise and process noise, as shown in thefollowing equation:

x(k+1)=Ax(k)+Bu(k)+w(k)

y(k)=Cx(k)+Du(k)+v(k)  (Equation 3)

where w(k) is the process noise and v(k) is the measurement noise.

State/disturbance estimator 520 can use the modified system model ofEquation 3 to estimate the system state {circumflex over (x)}(k) asshown in the following equation:

x(k+1|k)=A{circumflex over (x)}(k|k−1)+Bu(k)+K(y(k)−ŷ(k|k−1))

ŷ(k|k−1)=C{circumflex over (x)}(k|k−1)+Du(k)  (Equation 4)

where {circumflex over (x)}(k+1|k) is the estimated/predicted state attime step k+1 given the measurement at time step k, ŷ(k|k−1) is thepredicted output at time step k given the measurement at time step k−1,and K is the estimator gain. Equation 4 describes the prediction step ofa Kalman filter. The Kalman filter may be the optimal linear filterunder some ideal assumptions (e.g., the system is linear, no plant-modelmismatch, the process and measurement noise are follow a white noisedistribution with known covariance matrices).

In some embodiments, the estimated state vector {circumflex over(x)}(k+1|k), the output vector ŷ(k|k−1), and the input vector u(k) aredefined as follows:

${{\hat{x}\left( {k + 1} \middle| k \right)} = \begin{bmatrix}{{\hat{T}}_{ia}\left( {k + 1} \middle| k \right)} \\{{\hat{T}}_{m}\left( {k + 1} \middle| k \right)} \\{\hat{I}\left( {k + 1} \middle| k \right)}\end{bmatrix}},\mspace{14mu} {{\hat{y}\left( k \middle| {k - 1} \right)} = \begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\overset{.}{\hat{Q}}}_{HVAC}\left( k \middle| {k - 1} \right)}\end{bmatrix}},{{u(k)} = \begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}}$

where {circumflex over (T)}_(ia) is an estimate of the zone airtemperature T_(ia), {circumflex over (T)}_(m) is an estimate of the zonemass temperature T_(m), Î is an estimate of the integrating disturbanceI, {circumflex over ({dot over (Q)})}_(HVAC) is an estimate of theheating or cooling load provided by HVAC equipment 308, T_(sp) is thetemperature setpoint, and T_(oa) is the outdoor air temperature.

The state estimation described by Equation 4 is effective in addressingstationary noise. However, unmeasured correlated disturbances are alsocommon in practice. Correlated disturbances can lead to biasedestimates/predictions. In some embodiments, state/disturbance estimator520 accounts for disturbances by adding a disturbance state to thesystem model, as shown in the following equation:

$\begin{matrix}{\begin{bmatrix}{\hat{x}\left( {k + 1} \middle| k \right)} \\{\hat{d}\left( {k + 1} \middle| k \right)}\end{bmatrix} = {{\begin{bmatrix}A & B_{d} \\0 & I\end{bmatrix}\begin{bmatrix}{\hat{x}\left( {k + 1} \middle| k \right)} \\{\hat{d}\left( {k + 1} \middle| k \right)}\end{bmatrix}} + {\quad{{{\begin{bmatrix}B \\0\end{bmatrix}{u(k)}} + {\begin{bmatrix}K_{x} \\K_{d}\end{bmatrix}\left( {{y(k)} - {\hat{y}\left( k \middle| {k - 1} \right)}} \right)\mspace{20mu} {\hat{y}\left( k \middle| {k - 1} \right)}}} = {{\begin{bmatrix}C & C_{d}\end{bmatrix}\begin{bmatrix}{\hat{x}\left( k \middle| {k - 1} \right)} \\{\hat{d}\left( k \middle| {k - 1} \right)}\end{bmatrix}} + {{Du}(k)}}}}}} & \left( {{Equation}\mspace{14mu} 5} \right)\end{matrix}$

where {circumflex over (d)}(k) is the unmeasured disturbance at timestep k. The effect of adding the disturbance state is similar to that ofthe integral action of a proportional-integral controller.State/disturbance estimator 520 may select values of B_(d) and C_(d) andmodify the estimation problem as shown in Equation 5. In someembodiments, state/disturbance estimator 520 selects values of B_(d) andC_(d) such that the system model is an integrating input and/or outputdisturbance model.

In some embodiments, the choice of disturbance models may have asignificant impact on the ability of model predictive controller 302 toadequately estimate the lumped thermal building mass temperature T_(m).Accurate estimation of the thermal mass temperature T_(m) enables modelpredictive controller 302 to effectively use the solid building mass tostore thermal energy when preheating or precooling building zone 310 inorder to time-shift the thermal energy load and minimize the cost ofoperating HVAC equipment 308.

System Identifier

System identifier 510 can be configured to perform a systemidentification process to identify the values of the parameters in thebuilding thermal mass model and the HVAC load model. The functionsperformed by system identifier 510 may be performed when operating inthe parameter identification mode. The parameters identified by systemidentifier 510 may define the values of the system matrices A, B, C, andD and the estimator gain K. System identifier 510 can perform the systemidentification process during the commissioning phase of the controlsystem and/or after the detection of a substantial change in the modelparameters that has resulted in significant closed-loop performancedeterioration. In some embodiments, the system identification processperformed by system identifier 510 is the same or similar to the systemidentification process described in U.S. Pat. No. 9,235,657 titled“System Identification and Model Development” and granted Jan. 12, 2016.The entire disclosure of U.S. Pat. No. 9,235,657 is incorporated byreference herein.

System identification (SI) is the art and science of buildingmathematical models of dynamic systems from observed input-output data.In general, a model for a system may be developed using three paradigms:white box, grey box or black box models. These three systemidentification modeling paradigms are summarized in FIG. 6.

For applications that require high fidelity models, white box models maybe derived via detailed first-principles modeling approaches. Thisapproach often requires a high degree of knowledge of the physicalsystem and may result in a high-order nonlinear dynamic model. While theresulting model may capture most of the relevant physical behavior asystem, there is usually a very high engineering investment required todevelop such a model. Given the complexities of the resulting models,often white box models cannot be used in real-time predictivecontrollers owing to real-time computational restrictions. Nevertheless,white box models may help to evaluate control approaches in simulationto identify and mitigate potential pitfalls of the approach.

In a black box modeling paradigm, little a priori knowledge about themodel structure is assumed except for standard assumptions likelinearity and stability of the model. Black box methods fit the bestmodel to the input-output data. While there is little to no assumed apriori knowledge about the model, which offers a deal of flexibility,there are a few drawbacks to the approach. Given that little about themodel structure is assumed, black box models often tend to have a largeamount of parameters, which requires a large amount of data to properlyselect the model order and appropriate input-output representation. Tocapture the necessary amount of input-output data needed for black boxmethods, an undesirable length for the system identification experimentmay be required.

Finally, grey box modeling assumes that there is a known model structureoften derived through simplified first-principals or semi-physicalmodeling. The model structure is parameterized and parameter estimationtechniques are used to estimate the parameters. Owing to theincorporation of more a priori knowledge into grey box models, theresulting number of parameters for these types of models are often lessthan black box models. As a result of fewer parameters as well as thefact that the model structure (given it is valid) already encodesimportant relationships, less input-output data is typically needed tofit a good model with respect to that needed to fit a good model forblack box methods.

In some embodiments, system identifier 510 uses a grey box model torepresent HVAC equipment 308 and building zone 310 (i.e., the controlledsystem). In a grey box model system identification approach, a modelstructure

is selected, which is subsequently parameterized with a parameter θ∈

⊂

^(d) where

is the set of admissible model parameter values. The resulting possiblemodels are given by the set:

M={

(θ)|θ∈

}  (Equation 6)

To fit a proper parameter vector to the controlled system, systemidentifier 510 can collect input-output data. The input-output data mayinclude samples of the input vector u(k) and corresponding samples ofthe output vector y(k), as shown in the following equation:

Z ^(N)=[y(1),u(1),y(2),u(2), . . . ,y(N),u(N)]  (Equation 7)

where N>0 is the number of samples collected. Parameter estimation isthe problem of selecting the proper value {circumflex over (θ)}_(N) ofthe parameter vector given the input-output data Z^(N) (i.e., it is themapping: Z^(N)→{circumflex over (θ)}_(N)∈

).

In some embodiments, system identifier 510 uses a prediction errormethod to select the values of the parameter vector {circumflex over(θ)}_(N). A prediction error method (PEM) refers to a particular familyof parameter estimation methods that is of interest in the context ofthe present disclosure. The prediction error method used by systemidentifier 510 may include fitting the parameter vector {circumflex over(θ)}_(N) by minimizing some function of the difference between thepredicted output and observed output. For example, let ŷ(k, θ) be thepredicted output at time step k given the past input-output sequenceZ^(k-1) and the model

(θ) for θ∈

. The prediction error at time step k is given by:

ε(k,θ):=y(k)−ŷ(k,θ)  (Equation 8)

In some embodiments, system identifier 510 filters the prediction errorsequence through a stable linear filter and defines the followingprediction performance metric:

$\begin{matrix}{{V_{N}\left( {\theta,Z^{N}} \right)} = {\frac{1}{N}{\sum_{k = 1}^{N}{\left( {ɛ_{F}\left( {k,\theta} \right)} \right)}}}} & \left( {{Equation}\mspace{14mu} 9} \right)\end{matrix}$

where

(⋅) is the cost function (e.g., a positive definite function) andε_(F)(k, θ) is the filtered prediction error. The parameter estimate{circumflex over (θ)}_(N) of a prediction error method is then given by:

{circumflex over (θ)}_(N)={circumflex over (θ)}_(N)(Z ^(N))=arg

V _(N)(θ,Z ^(N))  (Equation 10)

In various embodiments, system identifier 510 can use any parameterestimation method that corresponds to Equation 10. For example, systemidentifier 510 can define a quadratic cost function that is the squareof the prediction errors, as shown in the following equation:

(ε(k,θ))=

(y(k)−ŷ(k|k−1,θ))=∥y(k)−ŷ(k|k−1,θ)∥₂ ²  (Equation 11)

where ŷ(k|k−1, θ) denotes the one-step ahead prediction of the outputusing the model

(θ). When the prediction errors are independently and identicallydistributed random variables from a normal distribution (i.e., theprocess and measurement noise is Gaussian white noise) and the modelbeing identified is linear, the cost function of Equation 11 is optimalfrom a statistical point-of-view.

Still referring to FIG. 5, system identifier 510 is shown to include athermal mass storage model generator 512, a HVAC load model generator514, and an estimator gain identifier 516. Thermal mass storage modelgenerator 512 can be configured to estimate or identify the parametersof a thermal mass storage model that describes the relationship betweenindoor air temperature T_(ia), building mass temperature T_(m), outdoorair temperature T_(oa), internal load {dot over (Q)}_(other), andsensible HVAC load {dot over (Q)}_(HVAC). The indoor air temperatureT_(ia) represents the temperature of the air within building zone 310.The building mass temperature T_(m) represents the temperature of thesolid objects within building zone 310. The internal load {dot over(Q)}_(other) represents heat generated internally within building zone310 by building occupants, electronics, resistive loads, and othersources of heat. The sensible HVAC load {dot over (Q)}_(HVAC) representsthe amount of heating or cooling applied to building zone 310 by HVACequipment 308.

HVAC load model generator 514 can be configured to generate a HVAC loadmodel that describes the relationship between indoor air temperatureT_(ia), the temperature setpoint T_(sp), and the sensible HVAC load {dotover (Q)}_(HVAC). In some embodiments, the HVAC load model describes thedynamics of HVAC equipment 308 when providing heating or cooling tobuilding zone 310. For example, the HVAC load model may represent HVACequipment 308 as a proportional-integral (PI) control system thatprovides the sensible HVAC load {dot over (Q)}_(HVAC) as a function ofthe indoor air temperature T_(ia) and the temperature setpoint T_(sp).

Estimator gain identifier 516 can be configured to calculate thestate/disturbance estimation gain (i.e., the estimator gain K). Theestimator gain K may be referred to in other contexts as the disturbancematrix or Kalman filter gain. The operation of thermal mass storagemodel generator 512, HVAC load model generator 514, and estimator gainidentifier 516 are described in greater detail below.

Thermal Mass Storage Model Generator

Thermal mass storage model generator 512 can be configured to estimateor identify the parameters of a thermal mass storage model thatdescribes the relationship between indoor air temperature T_(ia),building mass temperature T_(m), outdoor air temperature T_(oa),internal load {dot over (Q)}_(other), and sensible HVAC load {dot over(Q)}_(HVAC). In some embodiments, thermal mass storage model generator512 performs a multi-stage system identification process to generate theparameters of the thermal mass storage model. The stages of the systemidentification process may include input-output data collection,generating a parameterized building thermal model, and estimating (i.e.,fitting) the parameters of the parameterized building thermal modelusing a parameter identification algorithm.

The first stage of the system identification process is input-outputdata collection. The type of input-output data collected to fit themodel parameters may have a significant impact in the usefulness andaccuracy of the resulting model. In some embodiments, thermal massstorage model generator 512 manipulates the temperature setpoint T_(sp)in an open-loop fashion as to sufficiently and persistently excite thesystem. For example, thermal mass storage model generator 512 canpersistently excite (PE) a signal {s(k): k=0, 1, 2, . . . } withspectrum ϕ(ω) if ϕ_(s)(ω)>0 for almost all ω where:

$\begin{matrix}{{\varphi_{s}(\omega)} = {\sum_{m = {- \infty}}^{\infty}\left( {\lim_{N\rightarrow\infty}{\frac{1}{N}{\sum_{k = 1}^{N}{{s(k)}{s\left( {k - m} \right)}}}}} \right)}} & \left( {{Equation}\mspace{14mu} 12} \right)\end{matrix}$

Input-output data collected during closed-loop operation may not providean information-rich dataset with respect to system identificationbecause the inputs tend to be correlated to themselves as well as to theunmeasured disturbances as a result of closing the loop with acontroller. As a result input-output data obtained via closed-loopoperation may not meet the PE condition without adding perturbations tothe inputs similar to that used in extremum seeking control.

Instead of closed-loop identification, thermal mass storage modelgenerator 512 can perform a system identification experiment is to meetthe PE condition. In some embodiments, thermal mass storage modelgenerator 512 performs the system identification experiment whilebuilding zone 310 is occupied. Accordingly, it may be desirable todesign the system identification experiment to be as short as possiblein order to prevent or minimize disruptions with building occupants. Insome embodiments, thermal mass storage model generator 512 makes smallchanges to the temperature setpoints T_(sp) (e.g., changes of 3-5° F.)in order to have minimal impact on the comfort of the buildingoccupants.

In some embodiments, thermal mass storage model generator 512 generatesthe excitation signals using a pseudorandom binary sequence (PRBS). FIG.7 shows a graph 700 of an example excitation temperature setpoint signal702 over a one day period. Thermal mass storage model generator 512 canbe configured to generate and provide temperature setpoint signal 702 tosmart thermostat 100 during the system identification experiment.Temperature setpoint signal 702 can be generated a priori and tailoredto meet the comfort constraints. The minimum and maximum temperaturevalues of temperature setpoint signal 702 may be different depending onwhether HVAC equipment 308 are operating in a heating dominated orcooling dominated mode. The values of temperature setpoint signal 702may define the values of temperature setpoint variable T_(sp) in theinput data set used to fit the model parameters.

For variable speed systems, excitation signals generated from a PRBS maynot be suitable in practice owing to concerns over saturation effectsleading to nonlinearities in the experimental dataset. This issue can bemonitored and analyzed via high-fidelity building simulations andsetpoint experiments on buildings. Modifications to the excitationsignal generation methodology can be evaluated first in simulation tounderstand the technical trade-offs and limitations. Experiments onbuildings can also be performed to help assess the use of other types ofexcitation signals.

Thermal mass storage model generator 512 can generate a set of inputdata by modulating the temperature setpoint T_(sp) and recording thevalues of the temperature setpoint T_(sp) and the outside airtemperature T_(oa) at each time step during the system identificationperiod. The set of input data may include a plurality of values of thetemperature setpoint T_(sp) and the outside air temperature T_(oa).Thermal mass storage model generator 512 can generate a set of outputdata by monitoring and recording zone conditions at each time step ofthe system identification period. The zone conditions may include thezone air temperature T_(ia) and the heating/cooling load {dot over(Q)}_(HVAC). The set of output data may include a plurality of values ofthe zone air temperature T_(ia) and the heating/cooling load {dot over(Q)}_(HVAC).

The second stage of the system identification process is generating aparameterized building thermal model. The parameterized building thermalmodel may describe the relationship between zone temperature T_(ia) andsensible HVAC load {dot over (Q)}_(HVAC). In some embodiments, theparameterized building thermal model is based on a thermal circuitrepresentation of building zone 310. An example a thermal circuitrepresenting the heat transfer characteristics of building zone 310 isshown in FIG. 8.

Referring now to FIG. 8, a thermal circuit 800 representing the heattransfer characteristics of building zone 310 is shown, according tosome embodiments. Thermal circuit 800 is a two thermal resistance, twothermal capacitance (2R2C) control-oriented thermal mass model. Thethermal energy balance within thermal circuit 800 is given by thefollowing linear differential equations:

$\begin{matrix}{{C_{ia}{\overset{.}{T}}_{ia}} = {{\frac{1}{R_{mi}}\left( {T_{m} - T_{ia}} \right)} + {\frac{1}{R_{oi}}\left( {T_{oa} - T_{ia}} \right)} + {\overset{.}{Q}}_{HVAC} + {\overset{.}{Q}}_{other}}} & \left( {{Equation}\mspace{14mu} 13} \right)\end{matrix}$

where T_(ia) represents the indoor air temperature, T_(oa) representsthe outdoor air temperature, T_(m) represents the thermal masstemperature (i.e., the average temperature of solid objects in buildingzone 310), C_(m) represents thermal capacitance of the thermal mass,C_(ia) represents the indoor air thermal capacitance, R_(mi) representsthe thermal resistance between the indoor air and the thermal mass,R_(oi) represents the thermal resistance between the indoor air and theoutdoor air, {dot over (Q)}_(HVAC) represents the sensible heat added tobuilding zone 310 by HVAC equipment 308 (or removed from building zone310 if cooling is provided), and {dot over (Q)}_(other) represents theheat load disturbance (e.g., internal heat generation within buildingzone 310 via solar radiation, occupancy, electrical equipment, etc.).

In some embodiments, heat transfer between adjacent zones is modeled aspart of the heat load disturbance {dot over (Q)}_(other) and thus, notexplicitly accounted for in the model. In other embodiments, thermalinteractions between adjacent zones can be accounted for by modeling thetemperature of adjacent zones as additional temperature nodes and addinga thermal resistor between each adjacent zone and building zone 310. Forbuildings that have significant thermal coupling between zones (e.g.,large stores, open floorplan buildings, etc.), it may be desirable tomodel the thermal interactions between adjacent zones.

The thermal dynamic model of Equation 13 describes the heat transferbetween the indoor air and the thermal mass, the heat transfer throughthe building walls, the sensible heat/cooling duty of HVAC equipment 308{dot over (Q)}_(HVAC), and the internal heat load {dot over (Q)}_(other)on the space generated by the occupants, electrical plug and lighting,and solar irradiation. Throughout this disclosure, the model of Equation13 is referred to as the thermal model. The thermal model has beendemonstrated to yield acceptable prediction accuracy while managing thetrade-off between model complexity (e.g., number of model parameters)and model prediction accuracy. In some embodiments, the internal heatload {dot over (Q)}_(other) cannot be measured and is therefore treatedas an unmeasured, time-varying disturbance.

HVAC Load Model Generator

HVAC load model generator 514 can be configured to generate a HVAC loadmodel that describes the relationship between indoor air temperatureT_(ia), the temperature setpoint T_(sp), and the sensible HVAC load {dotover (Q)}_(HVAC). In some embodiments, the HVAC load model describes thedynamics of HVAC equipment 308 when providing heating or cooling tobuilding zone 310. For example, the HVAC load model may represent HVACequipment 308 as a proportional-integral (PI) control system thatprovides the sensible HVAC load {dot over (Q)}_(HVAC) as a function ofthe indoor air temperature T_(ia) and the temperature setpoint T_(sp).

In some embodiments, the sensible load {dot over (Q)}_(HVAC) provided byHVAC equipment 308 is a function of the temperature setpoint T_(sp) aswell as other factors including the heat gained/lost through thebuilding walls, the internal heat load {dot over (Q)}_(other), and theheat transferred to/removed from the building thermal mass. To provideor remove heat to a particular zone, various actions may be performed byHVAC equipment 308. For example, for a chilled water cooled buildingwith a central air handling unit (AHU) serving several variable airvolume (VAV) terminal boxes, a chiller may generate chilled water. Thechilled water can be pumped to the AHU cooling coil and air can beforced over the cooling coils of the AHU to cool the air to a supply airtemperature setpoint (e.g., 55° F.). The cooled air can be deliveredthrough the supply air duct to terminal VAV boxes at the location ofbuilding zone 310. The flow rate of cooled air can be adjusted via adamper such that building zone 310 reaches and/or stays at thetemperature setpoint T_(sp).

In some embodiments, the dynamics of the entire process of providingcooling or heating to building zone 310 is non-negligible on thetime-scale of interest for utilizing the energy storage of a buildingmass. To account for the dynamics of the heating/cooling process, HVACload model generator 514 can model the dynamics between the temperaturesetpoint T_(sp) and the HVAC sensible load {dot over (Q)}_(HVAC). Anexample of a model that can be model that can be generated by HVAC loadmodel generator 514 is the following proportional-integral (PI)controller model:

{dot over (Q)} _(HVAC,j) =K _(p,j)ε_(sp) +K _(I,j)∫₀ ^(t)ε_(sp)(s)ds

ε_(sp) =T _(sp,j) −T _(ia)  (Equation 14)

where j∈{clg, hlg} is the index that is used to denote either heating orcooling mode. In some embodiments, different sets of parameters aregenerated for the heating and cooling modes. In other words, HVAC loadmodel generator 514 can generate a first set of HVAC load modelparameters for the cooling mode and a second set of HVAC load modelparameters for the heating mode. In some embodiments, the heating andcooling load is constrained to the following set: {dot over(Q)}_(HVAC,c1g)∈[−{dot over (Q)}_(c1g,cap), 0] and {dot over(Q)}_(HVAC,htg)∈[0, {dot over (Q)}_(htg,cap)] where {dot over(Q)}_(HVAC)<0 indicates that HVAC equipment 308 are providing coolingand {dot over (Q)}_(HVAC)>0 indicates that HVAC equipment 308 areproviding heating

The model of Equation 14 includes a proportional gain term (i.e.,K_(p,j)ε_(sp)) that accounts for the error ε_(sp) between the indoor airtemperature T_(ia) and the setpoint T_(sp,j). The model of Equation 14also includes an integrator term (i.e., K_(I,j)∫₀ ^(t)ε_(sp)(s)ds) toaccount for the time-varying nature of the HVAC load {dot over(Q)}_(HVAC,j) required to force the indoor temperature T_(ia) to itssetpoint T_(sp,j). The integrator value may be correlated with theinternal heat load disturbance {dot over (Q)}_(other) and thus, may playthe role of an integrating disturbance model. In some embodiments, themodel of Equation 14 does not represent a particular PI controller, butrather is a lumped dynamic model describing the relationship between theHVAC load {dot over (Q)}_(HVAC,j) and the temperature setpoint T_(sp,j).In some embodiments, the HVAC load model defines the HVAC load {dot over(Q)}_(HVAC,j) as a function of the zone air temperature T_(ia), thetemperature setpoint T_(sp,j), the zone air humidity H_(ia), a zone airhumidity setpoint H_(sp,j), zone heat load {dot over (Q)}_(other),and/or any of a variety of factors that affect the amount of heating orcooling provided by HVAC equipment. It is contemplated that the HVACload model can include any combination of these or other variables invarious embodiments.

In some embodiments, HVAC load model generator 514 accounts for stagedHVAC equipment in the HVAC load model. For staged equipment, thesensible heating and cooling load {dot over (Q)}_(HVAC,j) can be modeledto take values in the countable set {dot over (Q)}_(HVAC,j)={{dot over(Q)}_(s1), {dot over (Q)}_(s2), . . . , {dot over (Q)}_(sn)}, where {dotover (Q)}_(si) i∈{1, . . . , n} is the HVAC load with the first throughith stage on. The staging down and staging up dynamics can be assumed tobe negligible. Alternatively, the staging down and staging up dynamicscan be captured using a low-order linear model.

The model of Equation 14 allows for continuous values of {dot over(Q)}_(HVAC,j) and may not be directly applicable to staged HVACequipment. However, HVAC load model generator 514 can compute thetime-averaged HVAC load {dot over (Q)}_(HVAC,j) over a time period,which gives a continuous version of the HVAC load {dot over(Q)}_(HVAC,j) needed to meet the temperature setpoint T_(sp,j). Forexample, HVAC load model generator 514 can apply a filtering techniqueto the equipment stage on/off signal to compute a time-averaged versionof the HVAC load {dot over (Q)}_(HVAC,j). The filtered version of theHVAC load {dot over (Q)}_(HVAC,j) can then be used as an input to themodel of Equation 14 to describe the relationship between the HVAC load{dot over (Q)}_(HVAC,j) and the temperature setpoint T_(sp,j).Throughout this disclosure, the variable {dot over (Q)}_(HVAC,j) is usedto represent both the continuous HVAC load when HVAC equipment 308 havea variable speed and to represent the filtered HVAC load when HVACequipment 308 are staged.

Incorporating the thermal and the HVAC load models together and writingthe system of equations as a linear system of differential equationsgives the following state space representation:

$\begin{bmatrix}{\overset{.}{T}}_{ia} \\{\overset{.}{T}}_{m} \\\overset{.}{I}\end{bmatrix} = {\begin{bmatrix}{\frac{1}{C_{ia}}\left( {K_{p,j} - \frac{1}{R_{mi}} - \frac{1}{R_{oi}}} \right)} & \frac{1}{C_{ia}R_{mi}} & {- \frac{K_{p,j}K_{I,j}}{C_{ia}}} \\\frac{1}{C_{m}R_{mi}} & {- \frac{1}{C_{m}R_{mi}}} & 0 \\{- 1} & 0 & 0\end{bmatrix}{\quad{{\begin{bmatrix}T_{ia} \\T_{m} \\I\end{bmatrix} + {\begin{bmatrix}{- \frac{K_{p,j}}{C_{ia}}} & \frac{1}{C_{ia}R_{oi}} \\0 & 0 \\1 & 0\end{bmatrix}\begin{bmatrix}T_{{sp},j} \\T_{oa}\end{bmatrix}} + {\begin{bmatrix}\frac{1}{C_{ia}} \\0 \\0\end{bmatrix}{{\overset{.}{Q}}_{other}\mspace{20mu}\begin{bmatrix}T_{ia} \\{\overset{.}{Q}}_{{HVAC},j}\end{bmatrix}}}} = {{\begin{bmatrix}1 & 0 & 0 \\{- K_{p,j}} & 0 & K_{I,j}\end{bmatrix}\begin{bmatrix}T_{ia} \\T_{m} \\I\end{bmatrix}} + {\begin{bmatrix}0 & 0 \\K_{p,j} & 0\end{bmatrix}\begin{bmatrix}T_{{sp},j} \\T_{oa}\end{bmatrix}}}}}}$

or more compactly:

{dot over (x)}=A _(c)(θ)x+B _(c)(θ)u+B _(d) d

y=C _(c)(θ)x+D _(c)(θ)u  (Equation 15)

where x^(T)=[T_(ia), T_(m), l], u^(T)=[T_(sp,j),T_(oa)], d={dot over(Q)}_(other)/C_(ia), y^(T)=[T_(ia), {dot over (Q)}_(HVAC,j)], θ is aparameter vector containing all non-zero entries of the system matrices,and

${{A_{c}(\theta)} = \begin{bmatrix}{- {\theta_{5}\left( {\theta_{1} + \theta_{2} + \theta_{4}} \right)}} & {\theta_{5}\theta_{2}} & {\theta_{3}\theta_{4}\theta_{5}} \\{\theta_{6}\theta_{2}} & {{- \theta_{6}}\theta_{2}} & 0 \\{- 1} & 0 & 0\end{bmatrix}},{{B_{c}(\theta)} = \begin{bmatrix}{\theta_{1}\theta_{5}} & {\theta_{4}\theta_{5}} \\0 & 0 \\1 & 0\end{bmatrix}},{B_{d} = \begin{bmatrix}\theta_{5} \\0 \\0\end{bmatrix}},{{C_{c}(\theta)} = \begin{bmatrix}1 & 0 & 0 \\{- \theta_{4}} & 0 & {\theta_{3}\theta_{4}}\end{bmatrix}},{{C_{d}(\theta)} = \begin{bmatrix}0 \\0\end{bmatrix}},{{D_{c}(\theta)} = \begin{bmatrix}0 & 0 \\\theta_{4} & 0\end{bmatrix}}$

The third stage of the system identification process includes fittingthe parameters θ of the parameterized building thermal model (i.e., themodel of Equation 15) to the input-output data using a parameteridentification algorithm. FIG. 9 is a flowchart of a process 900 whichcan be performed by system identifier 510 to fit the parameters θ of theparameterized building thermal model. If HVAC equipment 308 is staged,system identifier 510 can filter input-output data to remove the highfrequency dither in the indoor air temperature T_(ia) and can compute atime-averaged version of the HVAC load {dot over (Q)}_(HVAC) from thediscrete HVAC staging trajectory (step 902). System identifier 510 canthen fit the thermal model parameters and HVAC load model parameters θto the input-output data (or filtered input-output data) obtained duringthe first stage of the system identification process (step 904). Systemidentifier 510 can augment the resulting state-space model with anotherintegrating disturbance model and can estimate the Kalman filter gainfor the resulting model (step 906). Using data obtained in a secondaryexperiment or under normal operation, system identifier 510 can validatethe model through the use of statistics that capture the multi-stepprediction accuracy of the resulting model (step 908). Each of thesesteps is described in detail below.

Step 902 may include filtering the input-output data to remove a highfrequency dither in the indoor air temperature T_(ia) and computing atime-averaged version of the HVAC load {dot over (Q)}_(HVAC) from thediscrete HVAC staging trajectory. When HVAC equipment 308 are staged,the indoor air temperature trajectory may have a high frequency dithercaused by HVAC equipment 308 staging up and down. For example, whenoperating in a cooling mode, HVAC equipment 308 can be staged up (e.g.,by activating discrete chillers) to drive the indoor air temperatureT_(ia) to the temperature setpoint T_(sp) minus a deadband. When theindoor air temperature T_(ia) reaches the temperature setpoint T_(sp)minus the deadband, HVAC stages can be deactivated, which results in anincrease in the temperature T_(ia). The indoor air temperature T_(ia)may increase until it reaches the temperature setpoint T_(sp) plus thedeadband, at which time the HVAC stages can be activated again. Theobserved effect may be a high frequency dither or oscillation of thetemperature T_(ia) around the setpoint T_(sp) (shown in FIG. 10). Also,the HVAC load trajectory may take discrete values as stages are turnedon and off. The duty cycle of the on/off trajectory of the stages maydepend on the heat transfer to/from building zone 310, the error betweenthe setpoint and indoor air temperature, and the setpoint. Thus, a goodindication of the HVAC load {dot over (Q)}_(HVAC) on a continuous scaleis some time-average of the HVAC load trajectory or stage on/offtrajectories.

To remove the high frequency dither and to perform the time-averaging onthe HVAC load trajectory, system identifier 510 can filter theinput-output data obtained from the SI experiment. In some embodiments,system identifier 510 uses a first-order Savitzky-Golay filter (SGF) tofilter the input-output data. The SGF involves fitting a polynomial to aset of data samples and evaluating the resulting polynomial at a singlepoint within the approximation interval. The SGF is equivalent toperforming discrete convolution with a fixed impulse response. However,the SGF is a non-casual filter. To address this issue, system identifier510 can use the nearest past filtered measurements in the systemidentification process. For example, let Δt_(filter) be the filteringwindow of the SGF. For the filter to compute its smoothed estimate ofthe data at a given time step t, it can use data from t−Δt_(filter)/2 tot+Δt_(filter)/2. To provide the filtered measurements (indoor airtemperature T_(ia) and HVAC load {dot over (Q)}_(HVAC)) at a sample timet_(k), the filtered measurement from t_(k)−Δt_(filter)/2 can be used.

FIG. 10 is a pair of graphs 1000 and 1050 illustrating an example of theindoor air temperature T_(ia), trajectory and the HVAC load {dot over(Q)}_(HVAC) trajectory before filtering with the SGF. Line 1002represents the outdoor air temperature T_(oa) and line 1004 representsthe unfiltered indoor air temperature T_(ia). In steady-state, theindoor air temperature T_(ia), may vary around the temperature setpointwithin a region that is plus and minus the control deadband 1106. Line1052 represents the unfiltered sensible load {dot over (Q)}_(HVAC) forstaged HVAC equipment.

FIG. 11 is a pair of graphs 1100 and 1150 illustrating an example of theindoor air temperature T_(ia), trajectory and the HVAC load {dot over(Q)}_(HVAC) trajectory after filtering with the SGF. Line 1102represents the outdoor air temperature T_(oa) and line 1004 representsthe filtered indoor air temperature T_(ia). In some embodiments, thefiltering removes the variation of temperature about the setpoint tosmooth the temperature trajectory. Line 1152 represents the filteredsensible load {dot over (Q)}_(HVAC) for staged HVAC equipment and may bea time-averaged version of the unfiltered sensible load 1052 {dot over(Q)}_(HVAC) shown in FIG. 10.

Referring again to FIG. 9, step 904 may include fitting the thermalmodel parameters and HVAC load model parameters to the input-output data(or filtered input-output data) obtained during the first stage of thesystem identification process. In step 904, the time-varying heat loaddisturbance d may be omitted. System identifier 510 can identify theparameters for the model given by:

$\begin{matrix}{{{\overset{.}{x}(t)} = {{\begin{bmatrix}\theta_{1} & \theta_{2} & \theta_{3} \\\theta_{4} & {- \theta_{4}} & 0 \\{- 1} & 0 & 0\end{bmatrix}{x(t)}} + {\begin{bmatrix}{\theta_{5}\theta_{7}} & {\theta_{5}\theta_{6}} \\0 & 0 \\1 & 0\end{bmatrix}{u(t)}}}}{{y(t)} = {{\begin{bmatrix}1 & 0 & 0 \\{- \theta_{7}} & 0 & \theta_{8}\end{bmatrix}{x(t)}} + {\begin{bmatrix}0 & 0 \\\theta_{7} & 0\end{bmatrix}{u(t)}}}}} & \left( {{Equation}\mspace{14mu} 16} \right)\end{matrix}$

where the effect of the time-varying disturbance d is accounted for viathe integrator model.

Using the prediction error method (PEM) described previously, systemidentifier 510 can obtain an estimate of the model parameters θ. In someembodiments, system identifier 510 estimates the model parameters θusing the functions provided by the Matlab System Identificationtoolbox. For example, the identification can be performed using thefunction greyest with the initial state option set to estimate, thedisturbance model option set to none, and the estimation focus option isset to simulation. Other options may not be specifically set andtherefore default values can be used.

The initial condition of the model of Equation 16 may be unknown, butcan be from the data. The disturbance component (i.e., the Kalman gainestimate K) may not be included in the model of Equation 16. Theestimation focus can be set to prediction so that the algorithmautomatically computes the weighting function of the prediction errors.The stability estimation focus option may perform the same weighting asthe prediction option, but also enforces model stability. Under theparameterization of the model of Equation 16, the values of theresistance, capacitance, and PI parameters may be computed from theparameter vector θ if necessary. Given that the resulting PEMoptimization problem is a nonlinear, non-convex problem, the initialguess on the parameter values may have a significant impact on theparameter values obtained, and multiple executions of the PEM solverprovided a different initial guess each time may converge to differentlocal minima.

To perform step 904 autonomously, model verification can be performedbefore process 900 continues to 906. If any of the verification stepsare not satisfied, the PEM solver can be reinitialized with a differentrandom initial guess. The model verification steps that can be performedin step 904 may include ensuring stability by checking the eigenvaluesof the obtained A. In some embodiments, the model verification stepsinclude ensuring that the identified parameters are greater than zero.This is a condition arising from the physical meaningfulness of themodel. The model verification steps may include ensuring that thethermal capacitance of the air C_(ia) is less than the thermalcapacitance C_(m) of the building mass. In some embodiments, the modelverification steps include verifying that the matrix A is wellconditioned in order to avoid prediction problems and model sensitivityto noise. If the condition number of A is less than 10500, theconditioning of A may be deemed to be acceptable.

Once the obtained model satisfies all the checks, process 900 mayproceed to step 906. It is possible that the conditions described aboveare not ever satisfied. However, this could be due to a number ofreasons that may be difficult to autonomously diagnose and potentially,be associated with a problem beyond the limitations of the SI procedure.To address this problem, a limit can be set for the number of timesthrough the steps above. If the parameter acceptability criteria is notsatisfied after the maximum number of iterations through step 904,system identifier 510 can report an error for a human operator toinvestigate and address.

Estimator Gain Identifier

As discussed above, the model of Equation 16 may only account for themain thermal dynamics of the system and may not identify an appropriateestimator gain. In step 906, the model of Equation 16 can be augmentedwith another integrating disturbance model and the estimator gain can beidentified. In some embodiments, step 906 is performed by estimator gainidentifier 516. Estimator gain identifier 516 can be configured toaugment the model of Equation 16 to produce the augmented model:

$\begin{bmatrix}{\overset{.}{\hat{x}}(t)} \\{\overset{.}{\hat{d}}(t)}\end{bmatrix} = {{\begin{bmatrix}A_{c} & B_{d} \\0 & 0\end{bmatrix}\begin{bmatrix}{\hat{x}(t)} \\{\hat{d}(t)}\end{bmatrix}} + {\begin{bmatrix}B_{c} \\0\end{bmatrix}{u(t)}} + {\underset{= {:{K{(\varphi)}}}}{\underset{}{\begin{bmatrix}{K_{x}(\varphi)} \\{K_{d}(\varphi)}\end{bmatrix}}}\left( {{y(t)} - {\hat{y}(t)}} \right)}}$${\hat{y}(t)} = {{\begin{bmatrix}C_{c} & C_{d}\end{bmatrix}\begin{bmatrix}{\hat{x}(t)} \\{\hat{d}(t)}\end{bmatrix}} + {D_{c}{u(t)}}}$

where the parameters of A_(c), B_(c), C_(c) and D_(c) are the same asthe parameters θ determined in step 904. The disturbance model can beselected such that B_(d)=1 and C_(d)=0, which is an input integratingdisturbance model. K(ϕ) is the identified estimator gain parameterizedwith the parameter vector ϕ as shown in the following equation:

$\begin{matrix}{{{K_{x}(\varphi)} = \begin{bmatrix}\varphi_{1} & \varphi_{2} \\0 & 0 \\\varphi_{3} & \varphi_{4}\end{bmatrix}},{{K_{d}(\varphi)} = \begin{bmatrix}\varphi_{5} & \varphi_{6}\end{bmatrix}}} & \left( {{Equation}\mspace{14mu} 17} \right)\end{matrix}$

where the elements corresponding to the mass temperature (i.e., themiddle row of the matrix K_(x)(ϕ)) are set to zero. If these elementsare allowed to take non-zero values, poor estimation performance mayresult. Accordingly, estimator gain identifier 516 may set the elementscorresponding to the mass temperature T_(m) to zero. The values for theparameters ϕ can be obtained using the function greyest on the augmentedmodel, using the same input-output data as that used in step 904.

Estimator gain identifier 516 can be configured to generate thefollowing augmented matrices:

${A_{aug} = \begin{bmatrix}A_{c} & B_{d} \\0 & 0\end{bmatrix}},{B_{aug} = \begin{bmatrix}B_{c} \\0\end{bmatrix}},{C_{aug} = \begin{bmatrix}C_{c} & C_{d}\end{bmatrix}},{D_{aug} = D_{c}}$

Similar to the verifications performed in step 904, estimator gainidentifier 516 can check for the stability of the observer systemA_(aug)−KC_(aug), and its conditioning number. In this case, theconditioning number of A_(aug)−KC_(aug) may be deemed acceptable if itis less than 70,000.

Step 908 may include validating the augmented model generated in step906. As discussed above, system identifier 510 may select the optimalmodel parameters that fit the input-output data with respect to anobjective function that depends on the one-step error prediction. Instep 904, the PEM is the one-step ahead prediction without estimation(i.e., open-loop estimation), whereas step 906 represents the one-stepahead closed-loop estimation problem. However, within the context ofMPC, multi-step open-loop predictions can be used to select the optimalinput trajectory over the prediction horizon. Thus, the model validationperformed in step 908 may ensure that the estimator is tunedappropriately and that the multi-step open-loop prediction performedwith the identified model provides sufficient accuracy.

In some embodiments, system identifier 510 collects another input-outputdata set for model validation. This data may be collected during normaloperation or from another SI experiment. Using this validation data set,system identifier 510 can generate one or more metrics that quantify themulti-step prediction accuracy. One example of a metric which can begenerated in step 908 is the Coefficient of Variation Weighted MeanAbsolute Prediction Error (CVWMAPE). CVWMAPE is an exponentiallyweighted average of the absolute prediction error at each time step thatquantities the prediction error over time. System identifier 510 cancalculate the CVWMAPE as follows:

$\begin{matrix}{{{{CVWMAPE}(k)} = \frac{\sum_{i = k}^{k + N_{h} - 1}{e^{{- i}/N_{h}}{{{y(i)} - {\hat{y}\left( i \middle| k \right)}}}}}{\sum_{i = k}^{k + N_{h} + 1}{e^{{- i}/N_{h}}{{y(i)}}}}},{k = 0},1,2,\ldots} & \left( {{Equation}\mspace{14mu} 18} \right)\end{matrix}$

where N_(h)∈

_(>0) is the prediction horizon, y(i) is the measured output at timestep i, and ŷ(i|k) is the predicted output with the identified modelgiven a measurement at time step k and the input sequence u(k), u(k+1),. . . , u(i−1). For notational simplicity the variable y is used todenote one of the two outputs (i.e., in this subsection, y is a scalar).The CVWMAPE can be computed for both outputs.

Another prediction error to consider is prediction error with respect tothe q-step ahead prediction (q∈

_(≥0)). System identifier 510 can calculate the Coefficient of VariationRoot Mean Squared Prediction Error (CVRMSPE) is to evaluate the q-stepahead prediction error. For example, system identifier 510 can calculatethe CVRMSPE for a range of values of q from zero-step ahead predictionup to N_(h)-step ahead prediction. Given a set of measured output values{y(0), . . . , y(M)} for M∈

_(≥0), the CVRMSPE is given by:

$\begin{matrix}{{{CVRMSPE}(q)} = \frac{{RMSPE}(q)}{\overset{\_}{y}(q)}} & \left( {{Equation}\mspace{14mu} 19} \right)\end{matrix}$

where

$\begin{bmatrix}{{\hat{T}}_{ia}\left( {k + 1} \middle| k \right)} \\{{\hat{T}}_{m}\left( {k + 1} \middle| k \right)} \\{\hat{I}\left( {k + 1} \middle| k \right)} \\{\hat{d}\left( {k + 1} \middle| k \right)}\end{bmatrix} = {{A_{aug}\begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\hat{T}}_{m}\left( k \middle| {k - 1} \right)} \\{\hat{I}\left( k \middle| {k - 1} \right)} \\{\hat{d}\left( k \middle| {k - 1} \right)}\end{bmatrix}} + {B_{aug}\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}} + {\begin{bmatrix}K_{x} \\K_{d}\end{bmatrix}\left( {\begin{bmatrix}{T_{ia}(k)} \\{{\overset{.}{Q}}_{HVAC}(k)}\end{bmatrix} - \begin{bmatrix}{T_{ia}\left( k \middle| {k - 1} \right)} \\{{\overset{.}{\hat{Q}}}_{HVAC}\left( k \middle| {k - 1} \right)}\end{bmatrix}} \right)}}$

for all q∈{0, . . . , N_(h)−1}. The CVRMSPE helps identify theprediction error over the duration of the optimization period.

Once the system identification process is complete, system identifier510 can provide the system matrices A, B, C, and D with identifiedparameters θ and the estimator gain K with identified parameters ϕ tostate disturbance estimator 520 and predictive optimizer 522. Statedisturbance estimator 520 can use the system matrices A, B, C, and D andthe estimator gain K to compute an estimate of the current state{circumflex over (x)}(k) and unmeasured disturbance {circumflex over(d)}(k) at each sample time.

Load/Rate Predictor

Referring again to FIG. 5, model predictive controller 302 is shown toinclude a load/rate predictor 518. Load/rate predictor 518 can beconfigured to predict the heat load disturbance {circumflex over ({dotover (Q)})}_(other)(k) and the utility rates {circumflex over(r)}_(elec)(k) for each time step k (e.g., k=1 . . . n) of anoptimization period. Load/rate predictor 518 can provide the load andrate predictions to predictive optimizer 522 for use in determining theoptimal temperature setpoints T_(sp).

In some embodiments, load/rate predictor 518 predicts the value of{circumflex over ({dot over (Q)})}_(other)(k) at each time step k usinga history of disturbance estimates {circumflex over (d)}_(hist)generated by state/disturbance estimator 520. As described above, thestate/disturbance estimation problem is given by:

$\begin{matrix}{{\begin{bmatrix}{T_{ia}\left( k \middle| {k - 1} \right)} \\{{\overset{.}{\hat{Q}}}_{HVAC}\left( k \middle| {k - 1} \right)}\end{bmatrix} = {{C_{avg}\begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\hat{T}}_{m}\left( k \middle| {k - 1} \right)} \\{\hat{I}\left( k \middle| {k - 1} \right)} \\{\hat{d}\left( k \middle| {k - 1} \right)}\end{bmatrix}} + {D\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}}}}\mspace{20mu} {where}\mspace{20mu} {{A_{aug} = \begin{bmatrix}A & B_{d} \\0 & I\end{bmatrix}},\mspace{20mu} {B_{aug} = \begin{bmatrix}B \\0\end{bmatrix}},\mspace{20mu} {C_{aug} = {\begin{bmatrix}C & 0\end{bmatrix}.}}}} & \left( {{Equation}\mspace{14mu} 22} \right)\end{matrix}$

State/disturbance estimator 520 can perform the state/disturbanceestimation for each time step k to generate a disturbance estimate{circumflex over (d)}(k) for each time step k. State/disturbanceestimator 520 can record a history of disturbance estimates {circumflexover (d)}_(hist):={{circumflex over (d)}(k|k−1)}_(k=1) ^(N) ^(hist) andprovide the history of disturbance estimates {circumflex over(d)}_(hist) to load/rate predictor 518.

Load/rate predictor 518 can use the history of disturbance estimates{circumflex over (d)}_(hist) to predict the value of {circumflex over({dot over (Q)})}_(other)(k) at each time step k. In some embodiments,the heat load disturbance {circumflex over ({dot over (Q)})}_(other)(k)is a function of the time of day, the day type (e.g., weekend, weekday,holiday), the weather forecasts, building occupancy, and/or otherfactors which can be used by load/rate predictor 518 to predict the heatload disturbance {circumflex over ({dot over (Q)})}_(other)(k). Forexample, load/rate predictor 518 can predict the heat load disturbanceusing the equation:

{circumflex over ({dot over (Q)})}_(other)(k)=f _(pred)(kΔ+t ₀,{circumflex over (d)} _(hist))  (Equation 23)

where Δ>0 is the sample period and t₀ is an initial time. In someembodiments, load/rate predictor 518 uses weather forecasts from aweather service 526 and/or load history data from historical data 528 topredict the heat load disturbance {circumflex over ({dot over(Q)})}_(other)(k) at each time step.

In some embodiments, load/rate predictor 518 uses a deterministic plusstochastic model trained from historical load data to predict the heatload disturbance {circumflex over ({dot over (Q)})}_(other)(k).Load/rate predictor 518 may use any of a variety of prediction methodsto predict the heat load disturbance {circumflex over ({dot over(Q)})}_(other)(k) (e.g., linear regression for the deterministic portionand an AR model for the stochastic portion). In some embodiments,load/rate predictor 518 makes load/rate predictions using the techniquesdescribed in U.S. patent application Ser. No. 14/717,593.

Load/rate predictor 518 is shown receiving utility rates from utilities524. Utility rates may indicate a cost or price per unit of a resource(e.g., electricity, natural gas, water, etc.) provided by utilities 524at each time step k in the prediction window. In some embodiments, theutility rates are time-variable rates. For example, the price ofelectricity may be higher at certain times of day or days of the week(e.g., during high demand periods) and lower at other times of day ordays of the week (e.g., during low demand periods). The utility ratesmay define various time periods and a cost per unit of a resource duringeach time period. Utility rates may be actual rates received fromutilities 524 or predicted utility rates estimated by load/ratepredictor 518.

In some embodiments, the utility rates include demand charges for one ormore resources provided by utilities 524. A demand charge may define aseparate cost imposed by utilities 524 based on the maximum usage of aparticular resource (e.g., maximum energy consumption) during a demandcharge period. The utility rates may define various demand chargeperiods and one or more demand charges associated with each demandcharge period. In some instances, demand charge periods may overlappartially or completely with each other and/or with the predictionwindow. Predictive optimizer 522 may be configured to account for demandcharges in the high level optimization process performed by predictiveoptimizer 522. Utilities 524 may be defined by time-variable (e.g.,hourly) prices, a maximum service level (e.g., a maximum rate ofconsumption allowed by the physical infrastructure or by contract) and,in the case of electricity, a demand charge or a charge for the peakrate of consumption within a certain period. Load/rate predictor 518 maystore the predicted heat load disturbance {circumflex over ({dot over(Q)})}_(other)(k) and the utility rates in memory 508 and/or provide thepredicted heat load disturbance {circumflex over ({dot over(Q)})}_(other)(k) and the utility rates to predictive optimizer 522.

Predictive Optimizer

Still referring to FIG. 5, model predictive controller 302 is shown toinclude a predictive optimizer 522. Predictive optimizer 522 canmanipulate the temperature setpoint T_(sp) when operating in theoperational mode to minimize the economic cost of operating HVACequipment 308 over the duration of the optimization period. To determinethe temperature setpoints T_(sp), predictive optimizer 522 can optimizean objective function (i.e., a cost function) that accounts for the costof operating HVAC equipment 308 over the duration of the optimizationperiod. The costs of operating HVAC equipment 308 can include, forexample, the costs of resources consumed by HVAC equipment 308 duringoperation (e.g., electricity, natural gas, water, etc.), demand chargesimposed by an electric utility, peak load contribution charges,equipment degradation/replacement costs, and/or other costs associatedwith the operation of HVAC equipment 308.

An example of an objective function which can be optimized by predictiveoptimizer 522 is shown in the following equation:

min_(T) _(sp) _(,∈) _(T) _(,δT) _(sp) Σ_(k=0) ^(N-1)(r _(elec)(k)

_(elec)(k)+p _(∈) _(T) ∈_(T)(k)+p _(δT) δT _(sp)(k))   (Equation 24)

where r_(elec)(k) is the price of electricity at time step k,

_(elec)(k) is the predicted electricity consumption of HVAC equipment308 at time step k, ∈_(T)(k) is a number of degrees by which the zonetemperature constraints are violated at time step k, p_(∈) _(T) is azone air temperature penalty coefficient applied to ∈_(T)(k), δT_(sp)(k)is a number of degrees by which the zone temperature setpoint T_(sp)changes between time step k−1 and time step k, and p_(δT) _(sp) is atemperature setpoint change penalty coefficient applied to δT_(sp)(k).Equation 24 accounts for the cost of electricity consumption(r_(elec)(k)

_(elec)(k)), penalizes violations of the indoor air temperature bounds(p_(∈) _(T) ∈_(T)(k)), and penalizes changes in the temperature setpoint(p_(δT)δT_(sp)(k)).

If additional resources (other than electricity) are consumed by HVACequipment 308, additional terms can be added to Equation 24 to representthe cost of each resource consumed by HVAC equipment 308. For example,if HVAC equipment 308 consume natural gas and water in addition toelectricity, Equation 24 can be updated to include the terms r_(gas)(k):

_(gas)(k) and r_(water)(k)

_(water)(k) in the summation, where r_(gas)(k) is the cost per unit ofnatural gas at time step k,

_(gas)(k) is the predicted natural gas consumption of HVAC equipment 308at time step k, r_(water)(k) is the cost per unit of water at time stepk, and

_(water)(k) is the predicted water consumption of HVAC equipment 308 attime step k.

Predictive optimizer 522 can be configured to estimate the amount ofeach resource consumed by HVAC equipment 308 (e.g., electricity, naturalgas, water, etc.) as a function of the sensible heating or cooling load{circumflex over ({dot over (Q)})}_(HVAC). In some embodiments, aconstant efficiency model is used to compute the resource consumption ofHVAC equipment 308 as a fixed multiple heating or cooling load{circumflex over ({dot over (Q)})}_(HVAC) (e.g.,

_(elec)(k)=η{circumflex over ({dot over (Q)})}_(HVAC)(k)). In otherembodiments, equipment models describing a relationship between resourceconsumption and heating/cooling production can be used to approximatethe resource consumption of the HVAC equipment 308. The equipment modelsmay account for variations in equipment efficiency as a function of loadand/or other variables such as outdoor weather conditions (e.g.,

_(elec)(k)=f({circumflex over ({dot over (Q)})}_(HVAC)(k), T_(oa)(k),η_(base))). Predictive optimizer 522 can convert the estimatedheating/cooling load {circumflex over ({dot over (Q)})}_(HVAC)(k) ateach time step to one or more resource consumption values (e.g.,

_(elec)(k),

_(gas)(k),

_(water)(k), etc.) for inclusion in the objective function.

In some embodiments, predictive optimizer 522 is configured to modifythe objective function to account for various other costs associatedwith operating HVAC equipment 308. For example, predictive optimizer 522can modify the objective function to account for one or more demandcharges, peak load contribution charges, equipment degradation costs,equipment purchase costs, revenue generated from participating inincentive-based demand response programs, economic load demand response,and/or any other factors that can contribute to the cost incurred byoperating HVAC equipment 308 or revenue generated by operating HVACequipment 308. Predictive optimizer 522 can add one or more additionalterms to the objective function to account for these or other factors.Several examples of such functionality are described in U.S. patentapplication Ser. No. 15/405,236 filed Jan. 12, 2017, U.S. patentapplication Ser. No. 15/405,234 filed Jan. 12, 2017, U.S. patentapplication Ser. No. 15/426,962 filed Feb. 7, 2017, U.S. patentapplication Ser. No. 15/473,496 filed Mar. 29, 2017, and U.S. patentapplication Ser. No. 15/616,616 filed Jun. 7, 2017. The entiredisclosure of each of these patent applications is incorporated byreference herein. In some embodiments, model predictive controller 302includes some or all of the functionality of the controllers and/orcontrol systems described in these patent applications.

Predictive optimizer 522 can be configured to automatically generate andimpose constraints on the optimization of the objective function. Theconstraints may be based on the system model provided by systemidentifier 510, the state estimates provided by state/disturbanceestimator 520, the heat load disturbance predictions provided byload/rate predictor 518, and constraints on the zone air temperatureT_(ia). For example, predictive optimizer 522 can generate and imposethe following constraints for each time step k∈{0, . . . , N−1}:

$\begin{matrix}{\begin{bmatrix}{T_{ia}\left( {k + 1} \right)} \\{T_{m}\left( {k + 1} \right)} \\{I\left( {k + 1} \right)}\end{bmatrix} = {{A\begin{bmatrix}{T_{ia}(k)} \\{T_{m}(k)} \\{I(k)}\end{bmatrix}} + {B\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}} + {B_{d}{{\overset{.}{Q}}_{other}(k)}}}} & \left( {{Equation}\mspace{14mu} 25} \right) \\{\mspace{79mu} {\left\lbrack {Q_{HVAC}(k)} \right\rbrack = {{C\begin{bmatrix}{T_{ia}(k)} \\{T_{m}(k)} \\{I(k)}\end{bmatrix}} + {D\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}}}}} & \left( {{Equation}\mspace{14mu} 26} \right) \\{\mspace{85mu} {{{\hat{}}_{elec}(k)} = {\eta {{\overset{.}{\hat{Q}}}_{HVAC}(k)}}}} & \left( {{Equation}\mspace{14mu} 27} \right) \\{\mspace{79mu} {{T_{ia}(k)} \leq {{T_{\max}(k)} + T_{deadband} + {\epsilon_{T}(k)}}}} & \left( {{Equation}\mspace{14mu} 28} \right) \\{\mspace{79mu} {{T_{ia}(k)} \geq {{T_{\min}(k)} - T_{deadband} - {\epsilon_{T}(k)}}}} & \left( {{Equation}\mspace{14mu} 29} \right) \\{\mspace{79mu} {{T_{sp}(k)} \in \left\lbrack {{T_{{sp},\min}(k)} \leq {T_{{sp},\max}(k)}} \right\rbrack}} & \left( {{Equation}\mspace{14mu} 30} \right) \\{\mspace{79mu} {{{T_{sp}\left( {k + 1} \right)} - {T_{sp}(k)}} \leq {\delta \; {T_{sp}(k)}}}} & \left( {{Equation}\mspace{14mu} 31} \right) \\{\mspace{79mu} {{{- {T_{sp}\left( {k + 1} \right)}} + {T_{sp}(k)}} \leq {\delta \; {T_{sp}(k)}}}} & \left( {{Equation}\mspace{14mu} 32} \right) \\{\mspace{79mu} {{\epsilon_{T}(k)} \geq 0}} & \left( {{Equation}\mspace{14mu} 33} \right) \\{\mspace{79mu} {{{T_{ia}(0)} = {\hat{T}}_{{ia},0}},{{T_{m}(0)} = {\hat{T}}_{m,0}},{{I(0)} = {\hat{I}}_{0}}}} & \left( {{Equation}\mspace{14mu} 34} \right)\end{matrix}$

Equation 25 is based on the thermal mass storage model provided bythermal mass storage model generator 512. The thermal mass storage modelpredicts the system states at each time step k+1 as a function of thesystem states x(k), inputs u(k), and heat load disturbance at theprevious time step k. Specifically, Equation 25 defines the relationshipbetween the system states T_(ia)(k+1), T_(m)(k+1), and I(k+1) at timestep k+1, the system states T_(ia)(k), T_(m)(k), and I(k) at time stepk, the controlled or measured inputs T_(sp)(k) and T_(oa)(k) at timestep k, and the estimated heat load disturbance {dot over(Q)}_(other)(k) at time step k. The parameters in the matrices A and Bcan be determined by system identifier 510, as previously described. Thevalue of {dot over (Q)}_(other)(k) at each time step k can be determinedby load/rate predictor 518.

Equation 26 is based on the HVAC load model provided by HVAC load modelgenerator 514. The HVAC load model defines the heating/cooling load {dotover (Q)}_(HVAC)(k) at each time step k as a function of the systemstates T_(ia)(k), T_(m)(k), and I(k) at time step k and the controlledor measured inputs T_(sp)(k) and T_(oa)(k) at time step k. Theparameters in the matrices C and D can be determined by systemidentifier 510, as previously described.

Equation 27 defines the relationship between the predicted electricconsumption

_(elec)(k) at time step k and the heating/cooling load {dot over(Q)}_(HVAC)(k) at time step k. If additional resources (other thanelectricity) are consumed by HVAC equipment 308, additional constraintscan be added to define the relationship between {dot over (Q)}_(HVAC)(k)and the amount of each resource consumed by HVAC equipment 308.

Equations 28 and 29 constrain the zone air temperature T_(ia)(k) at eachtime step k between the minimum temperature threshold T_(min)(k) and themaximum temperature threshold T_(max)(k) at that time step k, plus orminus the temperature deadband T_(deadband) and a temperature error∈_(T)(k). In some embodiments, the values of the minimum temperaturethreshold T_(min)(k) and the maximum temperature threshold T_(max)(k)can vary over the duration of the optimization period. For example, theminimum temperature threshold T_(min)(k) and the maximum temperaturethreshold T_(max)(k) may define a first zone temperature range (i.e.,T_(min)(k) to T_(max)(k)) during time step k and a second zonetemperature range (i.e., T_(min)(k+1) to T_(max)(k+1)) during time stepk+1. The values of T_(min) and T_(max) at each time step can vary as afunction of the time of day, day of the week, building occupancy, orother factors. Predictive optimizer 522 may be allowed to violate thetemperature constraints if necessary, but any deviation from the definedtemperature range (i.e., between T_(min)(k)−T_(deadband) andT_(max)(k)+T_(deadband)) may be penalized in the objective function byimposing a penalty p_(∈) _(T) on the temperature error ∈_(T)(k).Equation 33 requires the temperature error ∈_(T)(k) to be non-negativeat each time step k.

Equation 30 limits the temperature setpoint T_(sp)(k) at each time stepbetween a minimum temperature setpoint T_(sp,min)(k) and a maximumtemperature setpoint T_(sp,max)(k), whereas Equations 31-32 define thechanges θT_(sp)(k) in the temperature setpoint T_(sp)(k) betweenconsecutive time steps. The minimum temperature setpoint T_(sp,min)(k)and the maximum temperature setpoint T_(sp,max)(k) may be different fromthe minimum temperature threshold T_(min)(k) and the maximum temperaturethreshold T_(max)(k) used to constrain the zone air temperatureT_(ia)(k). In some embodiments, the values of the minimum temperaturesetpoint T_(sp,min) (k) and the maximum temperature setpointT_(sp,max)(k) can vary over the duration of the optimization period. Forexample, the minimum temperature setpoint T_(sp,min)(k) and the maximumtemperature setpoint T_(sp,max)(k) may define a first temperaturesetpoint range (i.e., T_(sp,min) (k) to T_(sp,max)(k)) during time stepk and a second temperature setpoint range (i.e., T_(sp,min)(k+1) toT_(sp,max)(k+1)) during time step k+1. The values of T_(sp,min) andT_(sp,max) at each time step can vary as a function of the time of day,day of the week, building occupancy, or other factors. Changes in thetemperature setpoint T_(sp) may be penalized in the objective functionby imposing a penalty p_(θT) _(sp) on the temperature change δT_(sp)(k).

Equation 34 sets the initial system states for the zone air temperatureT_(ia)(0), the zone mass temperature T_(m)(0), and the integrateddisturbance I(0) at the first time step k=0 of the optimization period.The initial system states T_(ia)(0), T_(m)(0), and I(0) may be set tothe initial state estimates {circumflex over (T)}_(ia,0), {circumflexover (T)}_(m,0), and Î₀ provided by state/disturbance estimator 520. Inother words, the predictive model shown in Equations 25-26 may beinitialized using the state estimates provided by state/disturbanceestimator 520.

Predictive optimizer 522 can be configured to optimize the objectivefunction (Equation 24) subject to the optimization constraints(Equations 25-34) to determine optimal values of the temperaturesetpoint T_(sp)(k) at each time step k of the optimization period.Predictive optimizer 522 can use any of a variety of optimizationtechniques to perform the optimization. For example, the optimizationproblem can be formulated as a linear program (Equations 24-34) andsolved using a linear optimization technique (e.g., a basis exchangealgorithm, an interior point algorithm, a cutting-plane algorithm, abranch and bound algorithm, etc.) using any of a variety solvers and/orprogramming languages. In some embodiments, predictive optimizer 522performs the optimization at the beginning of each time step k todetermine optimal temperature setpoints T_(sp)(k) for the next N timesteps.

Once the optimal temperature setpoints T_(sp) have been determined,predictive optimizer 522 can provide the optimal temperature setpointT_(sp)(k) for the current time step k to smart thermostat 100 and/orequipment controller 406. State/disturbance estimator 520 can use thetemperature setpoint T_(sp)(k) for the current time step k along withfeedback received from HVAC equipment 308 and building zone 310 duringtime step k (e.g., measurements of T_(ia)(k), T_(oa)(k), {dot over(Q)}_(HVAC)(k), etc.) to update the state/disturbance estimates the nexttime the optimization is performed.

Model Predictive Control Processes

Referring now to FIGS. 12-14, several flowcharts of model predictivecontrol processes 1200-1400 are shown, according to some embodiments.Processes 1200-1400 can be performed by one or more components of modelpredictive control systems 300-400, as described with reference to FIGS.3-11. For example, processes 1200-1400 can be performed by modelpredictive controller 302, smart thermostat 100, equipment controller406, and/or HVAC equipment 308 to control the temperature of buildingzone 310.

Referring specifically to FIG. 12, process 1200 is shown to includeoperating a model predictive controller in system identification mode toidentify parameters of a thermal mass storage model and a HVAC loadmodel for a building zone (step 1202). In some embodiments, step 1202 isperformed by system identifier 510. Step 1202 can include generating athermal mass storage model by modeling the heat transfer characteristicsof the building zone using a thermal circuit 800, as described withreference to FIG. 8. Step 1202 may also include generating a HVAC loadmodel by modeling the amount of heating/cooling {dot over (Q)}_(HVAC)provided by HVAC equipment 308 as a function of the temperature setpointT_(sp) and the zone air temperature T_(ia). The thermal mass storagemodel and the HVAC load model can be used to generate a system ofdifferential equations that define the relationship between zone airtemperature T_(ia), zone mass temperature T_(m), outdoor air temperatureT_(oa), the amount of heating or cooling {dot over (Q)}_(HVAC) providedby HVAC equipment 308, the heat load disturbance {dot over (Q)}_(other),and the temperature setpoint T_(sp).

The unknown parameters in the thermal mass storage model and the HVACload model (e.g., R_(oi), R_(mi), C_(m), C_(ia)) can be organized intoparameter matrices A, B, C, and D, as shown in the following equations:

$\begin{bmatrix}{\overset{.}{T}}_{ia} \\{\overset{.}{T}}_{m} \\\overset{.}{I}\end{bmatrix} = {\begin{bmatrix}{\frac{1}{C_{ia}}\left( {K_{p,j} - \frac{1}{R_{mi}} - \frac{1}{R_{oi}}} \right)} & \frac{1}{C_{ia}R_{m\; i}} & {- \frac{K_{p,j}K_{I,j}}{C_{ia}}} \\\frac{1}{C_{m}R_{mi}} & {- \frac{1}{C_{m}R_{mi}}} & 0 \\{- 1} & 0 & 0\end{bmatrix}{\quad{{\begin{bmatrix}T_{ia} \\T_{m} \\I\end{bmatrix} + {\begin{bmatrix}{- \frac{K_{p,j}}{C_{ia}}} & \frac{1}{C_{ia}R_{oi}} \\0 & 0 \\1 & 0\end{bmatrix}\begin{bmatrix}T_{{sp},j} \\T_{oa}\end{bmatrix}} + {\begin{bmatrix}\frac{1}{C_{ia}} \\0 \\0\end{bmatrix}{{\overset{.}{Q}}_{other}\mspace{20mu}\begin{bmatrix}T_{ia} \\{\overset{.}{Q}}_{{HVAC},j}\end{bmatrix}}}} = {{\begin{bmatrix}1 & 0 & 0 \\{- K_{p,j}} & 0 & K_{I,j}\end{bmatrix}\begin{bmatrix}T_{ia} \\T_{m} \\I\end{bmatrix}} + {\begin{bmatrix}0 & 0 \\K_{p,j} & 0\end{bmatrix}\begin{bmatrix}T_{{sp},j} \\T_{oa}\end{bmatrix}}}}}}$

or more compactly:

{dot over (x)}=A _(c)(θ)x+B _(c)(θ)u+B _(d) d

y=C _(c)(θ)x+D _(c)(θ)u

where x^(T)=[T_(ia), T_(m), I], u^(T)=[T_(sp,j),T_(oa)], d={dot over(Q)}_(other)/C_(ia), y^(T)=[T_(ia), {dot over (Q)}_(HVAC,j)], θ is aparameter vector containing all non-zero entries of the system matrices,and

${{A_{c}(\theta)} = \begin{bmatrix}{- {\theta_{5}\left( {\theta_{1} + \theta_{2} + \theta_{4}} \right)}} & {\theta_{5}\theta_{2}} & {\theta_{3}\theta_{4}\theta_{5}} \\{\theta_{6}\theta_{2}} & {{- \theta_{6}}\theta_{2}} & 0 \\{- 1} & 0 & 0\end{bmatrix}},{{B_{c}(\theta)} = \begin{bmatrix}{\theta_{1}\theta_{5}} & {\theta_{4}\theta_{5}} \\0 & 0 \\1 & 0\end{bmatrix}},{B_{d} = \begin{bmatrix}\theta_{5} \\0 \\0\end{bmatrix}},{{C_{c}(\theta)} = \begin{bmatrix}1 & 0 & 0 \\{- \theta_{4}} & 0 & {\theta_{3}\theta_{4}}\end{bmatrix}}$ ${{C_{d}(\theta)} = \begin{bmatrix}0 \\0\end{bmatrix}},{{D_{c}(\theta)} = \begin{bmatrix}0 & 0 \\\theta_{4} & 0\end{bmatrix}}$

The unknown parameters θ in the thermal mass storage model and the HVACload model (i.e., the parameters in matrices A, B, C, and D) can beidentified by fitting the parameters θ to a set of input-output data(e.g., sets of values T_(ia), T_(oa), T_(sp), and {dot over (Q)}_(HVAC))collected from the HVAC system when operating in the systemidentification mode. In some embodiments, the input-output data arecollected by modulating the temperature setpoint T_(sp) and observingthe corresponding values of T_(ia), T_(oa), and {dot over (Q)}_(HVAC)for a plurality of time steps during the system identification mode. Theprocess of fitting the model parameters θ to the set of input-outputdata can be accomplished by performing process 900, as described withreference to FIG. 9.

If HVAC equipment 308 are staged, step 1202 can include filteringinput-output data to remove the high frequency dither in the indoor airtemperature T_(ia) and can compute a time-averaged version of the HVACload {dot over (Q)}_(HVAC) from the discrete HVAC staging trajectory.The thermal model parameters and HVAC load model parameters θ can thenbe fit to the input-output data (or filtered input-output data). In someembodiments, step 1202 includes augmenting the resulting state-spacemodel with another integrating disturbance model and estimating theKalman filter gain for the resulting model. Using data obtained in asecondary experiment or under normal operation, step 1202 can includevalidating the model through the use of statistics that capture themulti-step prediction accuracy of the resulting model.

Still referring to FIG. 12, process 1200 is shown to include operatingthe model predictive controller in operational mode to determinetemperature setpoints T_(sp) for the building zone (step 1204). In someembodiments, step 1204 is performed by predictive optimizer 522, asdescribed with reference to FIG. 5. To determine the temperaturesetpoints T_(sp), step 1204 can include optimizing an objective function(i.e., a cost function) that accounts for the cost of operating HVACequipment 308 over the duration of the optimization period. The costs ofoperating HVAC equipment 308 can include, for example, the costs ofresources consumed by HVAC equipment 308 during operation (e.g.,electricity, natural gas, water, etc.), demand charges imposed by anelectric utility, peak load contribution charges, equipmentdegradation/replacement costs, and/or other costs associated with theoperation of HVAC equipment 308. The objective function can be optimizedsubject to constraints provided by the thermal mass storage model andthe HVAC load model, as well as constraints on the zone temperatureT_(ia). Step 1204 can be accomplished by performing process 1300,described in greater detail with reference to FIG. 13.

In some embodiments, process 1200 includes automatically switchingbetween the system identification mode and the operational mode based ona prediction error of the thermal mass storage model and/or the HVACload model. For example, process 1200 can include comparing the value ofthe zone temperature {circumflex over (T)}_(ia)(k|k−1) predicted by thethermal mass storage model to the actual measured value of the zonetemperature T_(ia)(k). If the difference between the predicted value{circumflex over (T)}_(ia)(k|k−1) and the actual value T_(ia)(k) exceedsan error threshold, process 1200 may automatically return to step 1202and repeat the system identification process. Similarly, process 1200can include comparing the value of the HVAC equipment load {circumflexover ({dot over (Q)})}_(HVAC)(k|k−1) predicted by the HVAC load model tothe actual value of the HVAC load {dot over (Q)}_(HVAC)(k). If thedifference between the predicted value {circumflex over ({dot over(Q)})}_(HVAC)(k|k−1) and the actual value {dot over (Q)}_(HVAC)(k)exceeds an error threshold, process 1200 may automatically return tostep 1202 and repeat the system identification process.

Process 1200 is shown to include providing the optimal temperaturesetpoints T_(sp) to an equipment controller (step 1206) and using theequipment controller to operate HVAC equipment to drive the temperatureof the building zone T_(ia) to the optimal temperature setpoints T_(sp)(step 1208). Step 1206 can include sending the optimal temperaturesetpoints T_(sp) from model predictive controller 302 to an equipmentcontroller 406. In some embodiments, both model predictive controller302 and equipment controller 406 are components of a smart thermostat100 (as shown in FIG. 4). In other embodiments, model predictivecontroller 302 is separate from smart thermostat 100 and configured tosend the optimal temperature setpoints to smart thermostat 100 via acommunications network 304 (as shown in FIG. 3). Smart thermostat 100and/or equipment controller 406 can use the optimal temperaturesetpoints T_(sp) to generate control signals for HVAC equipment 308which operate to drive the temperature of the building zone T_(ia) tothe optimal temperature setpoints T_(sp).

Referring now to FIG. 13, a flowchart of a process 1300 for generatingoptimal temperature setpoints for a building zone is shown, according tosome embodiments. In some embodiments, process 1300 is performed bymodel predictive controller 302, as described with reference to FIGS.5-11. Process 1300 can be performed to accomplish step 1204 of process1200.

Process 1300 is shown to include predicting a heat load disturbance {dotover (Q)}_(other) experienced by a building zone at each time step (k=0. . . N−1) of an optimization period (step 1302). In some embodiments,step 1302 is performed by load/rate predictor 518 and state/disturbancepredictor 520. Step 1302 can include generating a disturbance estimate{circumflex over (d)}(k) for each time step k in the optimizationperiod. In some embodiments, step 1302 includes estimating thedisturbance state {circumflex over (d)}(k) at each time step k in theoptimization period using the following state/disturbance model:

$\begin{bmatrix}{{\hat{T}}_{ia}\left( {k + 1} \middle| k \right)} \\{{\hat{T}}_{m}\left( {k + 1} \middle| k \right)} \\{\hat{I}\left( {k + 1} \middle| k \right)} \\{\hat{d}\left( {k + 1} \middle| k \right)}\end{bmatrix} = {{{A_{aug}\begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\hat{T}}_{m}\left( k \middle| {k - 1} \right)} \\{\hat{I}\left( k \middle| {k - 1} \right)} \\{\hat{d}\left( k \middle| {k - 1} \right)}\end{bmatrix}} + {B_{aug}\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}} + {\begin{bmatrix}K_{x} \\K_{d}\end{bmatrix}{\left( {\begin{bmatrix}{T_{ia}(k)} \\{{\overset{.}{Q}}_{HVAC}(k)}\end{bmatrix} - \begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\overset{.}{\hat{Q}}}_{HVAC}\left( k \middle| {k - 1} \right)}\end{bmatrix}} \right)\mspace{20mu}\begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\overset{.}{\hat{Q}}}_{HVAC}\left( k \middle| {k - 1} \right)}\end{bmatrix}}}} = {{C_{aug}\begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\hat{T}}_{m}\left( k \middle| {k - 1} \right)} \\{\hat{I}\left( k \middle| {k - 1} \right)} \\{\hat{d}\left( k \middle| {k - 1} \right)}\end{bmatrix}} + {D\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}}}}$$\mspace{20mu} {{{{where}\mspace{14mu} A_{aug}} = \begin{bmatrix}A & B_{d} \\0 & I\end{bmatrix}},{B_{aug} = \begin{bmatrix}B \\0\end{bmatrix}},{C_{aug} = {\begin{bmatrix}C & 0\end{bmatrix}.}}}$

Step 1302 can include recording a history of disturbance estimates{circumflex over (d)}_(hist):={{circumflex over (d)}(k|k−1)}_(k=1) ^(N)^(hist) predicted by the state/disturbance model.

Step 1302 can include using the history of disturbance estimates{circumflex over (d)}_(hist) to predict the value of {circumflex over({dot over (Q)})}_(other)(k) at each time step k. In some embodiments,the heat load disturbance {circumflex over ({dot over (Q)})}_(other)(k)is a function of the time of day, the day type (e.g., weekend, weekday,holiday), the weather forecasts, building occupancy, and/or otherfactors which can be used by load/rate predictor 518 to predict the heatload disturbance {circumflex over ({dot over (Q)})}_(other)(k). Forexample, step 1302 can include predicting the heat load disturbanceusing the equation:

{circumflex over ({dot over (Q)})}_(other)(k)=f _(pred)(kΔ+t ₀,{circumflex over (d)} _(hist))

where Δ>0 is the sample period and t₀ is an initial time. In someembodiments, step 1302 includes using weather forecasts from a weatherservice 526 and/or load history data from historical data 528 to predictthe heat load disturbance {circumflex over ({dot over (Q)})}_(other)(k)at each time step.

In some embodiments, step 1302 includes using a deterministic plusstochastic model trained from historical load data to predict the heatload disturbance {circumflex over ({dot over (Q)})}_(other)(k). Step1302 can include using any of a variety of prediction methods to predictthe heat load disturbance {circumflex over ({dot over (Q)})}_(other)(k)(e.g., linear regression for the deterministic portion and an AR modelfor the stochastic portion). In some embodiments, step 1302 includespredicting the heat load disturbance {circumflex over ({dot over(Q)})}_(other)(k) using the techniques described in U.S. patentapplication Ser. No. 14/717,593.

Still referring to FIG. 13, process 1300 is shown to include estimatinginitial system states {circumflex over (x)}(0) for a first time step(k=0) of the optimization period (step 1304). The initial system states{circumflex over (x)}(0) may include the estimated zone air temperature{circumflex over (T)}_(ia)(0), the estimated zone mass temperature{circumflex over (T)}_(m) (0), and the estimated integrating disturbanceÎ(0) at the beginning of the optimization period. In some embodiments,step 1304 is performed by state/disturbance estimator 520. The initialsystem states {circumflex over (x)}(0) can be estimated using same thestate/disturbance estimation model used in step 1302. For example, step1304 can include estimating the system states {circumflex over (x)}(k)at each time step k using the following equation:

{circumflex over (x)}(k+1|k)=A{circumflex over(x)}(k|k−1)+Bu(k)+K(y(k)−ŷ(k|k−1))

ŷ(k|k−1)=C{circumflex over (x)}(k|k−1)+Du(k)

where {circumflex over (x)}(k+1|k) is the estimated/predicted state attime step k+1 given the measurement at time step k, ŷ(k|k−1) is thepredicted output at time step k given the measurement at time step k−1,and K is the estimator gain.

In some embodiments, the estimated state vector {circumflex over(x)}(k+1|k), the output vector ŷ(k|k−1), and the input vector u(k) aredefined as follows:

${{\hat{x}\left( {k + 1} \middle| k \right)} = \begin{bmatrix}{{\hat{T}}_{ia}\left( {k + 1} \middle| k \right)} \\{{\hat{T}}_{m}\left( {k + 1} \middle| k \right)} \\{\hat{I}\left( {k + 1} \middle| k \right)}\end{bmatrix}},{{\hat{y}\left( k \middle| {k - 1} \right)} = \begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\overset{.}{\hat{Q}}}_{HVAC}\left( k \middle| {k - 1} \right)}\end{bmatrix}},{{u(k)} = \begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}}$

where {circumflex over (T)}_(ia) is an estimate of the zone airtemperature T_(ia), {circumflex over (T)}_(m) is an estimate of the zonemass temperature T_(m), Î is an estimate of the integrating disturbanceI, {circumflex over ({dot over (Q)})}_(HVAC) is an estimate of theheating or cooling load provided by HVAC equipment 308, T_(sp) is thetemperature setpoint, and T_(oa) is the outdoor air temperature. Theestimated system states {circumflex over (x)}(0) at the first time stepk=0 can be used as the initial system states.

Still referring to FIG. 13, process 1300 is shown to include using athermal mass storage model and a heat load model to predict thetemperature of the building zone {circumflex over (T)}_(ia) and the HVACequipment load {circumflex over ({dot over (Q)})}_(HVAC) at each timestep as a function of the estimated system states {circumflex over (x)},temperature setpoint trajectory T_(sp), and estimated heat loaddisturbance {circumflex over ({dot over (Q)})}_(other) (step 1306). Thethermal mass storage model may be defined as follows:

$\begin{bmatrix}{{\hat{T}}_{ia}\left( {k + 1} \middle| k \right)} \\{{\hat{T}}_{m}\left( {k + 1} \middle| k \right)} \\{\hat{I}\left( {k + 1} \middle| k \right)}\end{bmatrix} = {{A\begin{bmatrix}{{\hat{T}}_{ia}(k)} \\{{\hat{T}}_{m}(k)} \\{\hat{I}(k)}\end{bmatrix}} + {B\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}} + {B_{d}{{\hat{\overset{.}{Q}}}_{other}(k)}}}$

where {circumflex over (T)}_(ia)(k+1) is the predicted temperature ofthe building zone at time step k+1, {circumflex over (T)}_(m)(k+1) isthe predicted temperature of the building mass at time step k+1, Î(k+1)is the predicted value of the integrating disturbance at time step k+1,{circumflex over (T)}_(ia)(k) is the predicted temperature of thebuilding zone at time step k, {circumflex over (T)}_(m)(k) is thepredicted temperature of the building mass at time step k, Î(k) is thepredicted value of the integrating disturbance at time step k, T_(sp)(k)is the temperature setpoint at time step k, T_(oa)(k) is the outdoor airtemperature (measured) at time step k, and {dot over ({circumflex over(Q)})}_(other)(k) is the estimated heat load disturbance at time step k.

The HVAC load model may be defined as follows:

$\left\lbrack {{\hat{\overset{.}{Q}}}_{HVAC}(k)} \right\rbrack = {{C\begin{bmatrix}{{\hat{T}}_{ia}(k)} \\{{\hat{T}}_{m}(k)} \\{\hat{I}(k)}\end{bmatrix}} + {D\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}}}$

where {dot over ({circumflex over (Q)})}_(HVAC)(k) is the predicted HVACequipment load at time step k and {circumflex over (T)}_(ia)(k),{circumflex over (T)}_(m)(k), Î(k), T_(sp)(k), and T_(oa)(k) are thesame as the corresponding variables in the thermal mass storage model.

Process 1300 is shown to include optimizing an economic cost functionsubject to constraints on the predicted temperature of the building zone{circumflex over (T)}_(ia) to determine optimal temperature setpointsT_(sp) for the building zone at each time step of the optimizationperiod (step 1308). An example of an objective function which can beoptimized in step 1308 is shown in the following equation:

$\min\limits_{T_{sp},\epsilon_{T},{\delta \; T_{sp}}}{\sum\limits_{k = 0}^{N - 1}\left( {{{r_{elec}(k)}{{\hat{}}_{elec}(k)}} + {p_{\epsilon_{T}}{\epsilon_{T}(k)}} + {p_{\delta \; T}\delta \; {T_{sp}(k)}}} \right)}$

where r_(elec)(k) is the price of electricity at time step k,

_(elec)(k) is the predicted electricity consumption of HVAC equipment308 at time step k, ∈_(T)(k) is a number of degrees by which the zonetemperature constraints are violated at time step k, p_(∈) _(T) is azone air temperature penalty coefficient applied to ∈_(T)(k), δT_(sp)(k)is a number of degrees by which the zone temperature setpoint T_(sp)changes between time step k−1 and time step k, and p_(δT) _(sp) is atemperature setpoint change penalty coefficient applied to δT_(sp)(k).Equation 24 accounts for the cost of electricity consumption(r_(elec)(k)

_(elec)(k)), penalizes violations of the indoor air temperature bounds(p_(∈) _(T) ∈_(T)(k)), and penalizes changes in the temperature setpoint(p_(δT)δT_(sp)(k)).

If additional resources (other than electricity) are consumed by HVACequipment 308, additional terms can be added to the objective functionto represent the cost of each resource consumed by HVAC equipment 308.For example, if HVAC equipment 308 consume natural gas and water inaddition to electricity, the objective function can be updated toinclude the terms r_(gas)(k)

_(gas)(k) and r_(water)(k)

_(water)(k) in the summation, where r_(gas)(k) is the cost per unit ofnatural gas at time step k,

_(gas)(k) is the predicted natural gas consumption of HVAC equipment 308at time step k, r_(water)(k) is the cost per unit of water at time stepk, and

_(water)(k) is the predicted water consumption of HVAC equipment 308 attime step k.

In some embodiments, step 1308 includes estimating the amount of eachresource consumed by HVAC equipment 308 (e.g., electricity, natural gas,water, etc.) as a function of the sensible heating or cooling load{circumflex over ({dot over (Q)})}_(HVAC). In some embodiments, aconstant efficiency model is used to compute the resource consumption ofHVAC equipment 308 as a fixed multiple heating or cooling load{circumflex over ({dot over (Q)})}_(HVAC) (e.g.,

_(elec)(k)=η{circumflex over ({dot over (Q)})}_(HVAC)(k)). In otherembodiments, equipment models describing a relationship between resourceconsumption and heating/cooling production can be used to approximatethe resource consumption of the HVAC equipment 308. The equipment modelsmay account for variations in equipment efficiency as a function of loadand/or other variables such as outdoor weather conditions (e.g.,

_(elec)(k)=f({circumflex over ({dot over (Q)})}_(HVAC)(k), T_(oa)(k),η_(base))). Step 1308 can include converting the estimatedheating/cooling load {circumflex over ({dot over (Q)})}_(HVAC)(k) ateach time step to one or more resource consumption values (e.g.,

_(elec) (k),

_(gas)(k),

_(water)(k), etc.) for inclusion in the objective function.

Step 1308 can include optimizing the objective function to determineoptimal values of the temperature setpoint T_(sp)(k) for each time stepin the optimization period. The objective function can be optimizedsubject to a set of constraints (e.g., Equations 25-34). The constraintsmay include the thermal mass storage model and the HVAC load model,which define the relationship between {circumflex over ({dot over(Q)})}_(HVAC), the temperature setpoints T_(sp), and the zone airtemperature T_(ia) at each time step. The constraints may also includeconstraints on the zone air temperature T_(ia)(k) and constraints thatdefine the penalty terms in the objective function.

In some embodiments, step 1308 includes performing the optimization atthe beginning of each time step k to determine optimal temperaturesetpoints T_(sp)(k) for the next N time steps. Once the optimaltemperature setpoints T_(sp) have been determined, the optimaltemperature setpoint T_(sp)(k) for the current time step k can beprovided to smart thermostat 100 and/or equipment controller 406 and theoptimization period can be shifted forward in time by one time step. Thetemperature setpoint T_(sp)(k) for the current time step k along withfeedback received from HVAC equipment 308 and building zone 310 duringtime step k (e.g., measurements of T_(ia)(k), T_(oa)(k), {dot over(Q)}_(HVAC)(k), etc.) can be used to update the state/disturbanceestimates the next time the optimization is performed. Steps 1304-1308can be repeated at the beginning of each time step to determine theoptimal setpoint trajectory T_(sp) for the new (i.e., time shifted)optimization period.

Referring now to FIG. 14, a flowchart of a process 1400 for controllingthe temperature of a building zone using model predictive control isshown, according to some embodiments. Process 1400 can be performed byone or more components of model predictive control systems 300-400, asdescribed with reference to FIGS. 3-11. For example, process 1400 can beperformed by model predictive controller 302, smart thermostat 100,equipment controller 406, and/or HVAC equipment 308 to control thetemperature of building zone 310.

Process 1400 is shown to include identifying parameters of a thermalmass storage model and a HVAC load model for a building zone (step1402). In some embodiments, step 1202 is performed by system identifier510. Step 1402 can include generating a thermal mass storage model bymodeling the heat transfer characteristics of the building zone using athermal circuit 800, as described with reference to FIG. 8. Step 1402may also include generating a HVAC load model by modeling the amount ofheating/cooling {dot over (Q)}_(HVAC) provided by HVAC equipment 308 asa function of the temperature setpoint T_(sp) and the zone airtemperature T_(ia). The thermal mass storage model and the HVAC loadmodel can be used to generate a system of differential equations thatdefine the relationship between zone air temperature T_(ia), zone masstemperature T_(m), outdoor air temperature T_(oa), the amount of heatingor cooling {dot over (Q)}_(HVAC) provided by HVAC equipment 308, theheat load disturbance {dot over (Q)}_(other), and the temperaturesetpoint T_(sp).

The unknown parameters in the thermal mass storage model and the HVACload model (e.g., R_(oi), R_(mi), C_(m), C_(ia)) can be organized intoparameter matrices A, B, C, and D, as shown in the following equations:

$\begin{bmatrix}{\overset{.}{T}}_{ia} \\{\overset{.}{T}}_{m} \\\overset{.}{I}\end{bmatrix} = {\begin{bmatrix}{\frac{1}{C_{ia}}\left( {K_{p,j} - \frac{1}{R_{mi}} - \frac{1}{R_{oi}}} \right)} & \frac{1}{C_{ia}R_{m\; i}} & {- \frac{K_{p,j}K_{I,j}}{C_{ia}}} \\\frac{1}{C_{m}R_{mi}} & {- \frac{1}{C_{m}R_{mi}}} & 0 \\{- 1} & 0 & 0\end{bmatrix}{\quad{{\begin{bmatrix}T_{ia} \\T_{m} \\I\end{bmatrix} + {\begin{bmatrix}{- \frac{K_{p,j}}{C_{ia}}} & \frac{1}{C_{ia}R_{oi}} \\0 & 0 \\1 & 0\end{bmatrix}\begin{bmatrix}T_{{sp},j} \\T_{oa}\end{bmatrix}} + {\begin{bmatrix}\frac{1}{C_{ia}} \\0 \\0\end{bmatrix}{{\overset{.}{Q}}_{other}\mspace{20mu}\begin{bmatrix}T_{ia} \\{\overset{.}{Q}}_{{HVAC},j}\end{bmatrix}}}} = {{\begin{bmatrix}1 & 0 & 0 \\{- K_{p,j}} & 0 & K_{I,j}\end{bmatrix}\begin{bmatrix}T_{ia} \\T_{m} \\I\end{bmatrix}} + {\begin{bmatrix}0 & 0 \\K_{p,j} & 0\end{bmatrix}\begin{bmatrix}T_{{sp},j} \\T_{oa}\end{bmatrix}}}}}}$

or more compactly:

{dot over (x)}=A _(c)(θ)x+B _(c)(θ)u+B _(d) d

y=C _(c)(θ)x+D _(c)(θ)u

where x^(T)=[T_(ia), T_(m), I], u^(T)=[T_(sp,j),T_(oa)], d={dot over(Q)}_(other)/C_(ia), y^(T)=[T_(ia), {dot over (Q)}_(HVAC,j)], θ is aparameter vector containing all non-zero entries of the system matrices,and

${{A_{c}(\theta)} = \begin{bmatrix}{- {\theta_{5}\left( {\theta_{1} + \theta_{2} + \theta_{4}} \right)}} & {\theta_{5}\theta_{2}} & {\theta_{3}\theta_{4}\theta_{5}} \\{\theta_{6}\theta_{2}} & {{- \theta_{6}}\theta_{2}} & 0 \\{- 1} & 0 & 0\end{bmatrix}},{{B_{c}(\theta)} = \begin{bmatrix}{\theta_{1}\theta_{5}} & {\theta_{4}\theta_{5}} \\0 & 0 \\1 & 0\end{bmatrix}},{B_{d} = \begin{bmatrix}\theta_{5} \\0 \\0\end{bmatrix}},{{C_{c}(\theta)} = \begin{bmatrix}1 & 0 & 0 \\{- \theta_{4}} & 0 & {\theta_{3}\theta_{4}}\end{bmatrix}}$ ${{C_{d}(\theta)} = \begin{bmatrix}0 \\0\end{bmatrix}},{{D_{c}(\theta)} = \begin{bmatrix}0 & 0 \\\theta_{4} & 0\end{bmatrix}}$

The unknown parameters θ in the thermal mass storage model and the HVACload model (i.e., the parameters in matrices A, B, C, and D) can beidentified by fitting the parameters θ to a set of input-output data(e.g., sets of values T_(ia), T_(oa), T_(sp), and {dot over (Q)}_(HVAC))collected from the HVAC system when operating in the systemidentification mode. In some embodiments, the input-output data arecollected by modulating the temperature setpoint T_(sp) and observingthe corresponding values of T_(ia), T_(oa), and {dot over (Q)}_(HVAC)for a plurality of time steps during the system identification mode. Theprocess of fitting the model parameters θ to the set of input-outputdata can be accomplished by performing process 900, as described withreference to FIG. 9.

If HVAC equipment 308 are staged, step 1402 can include filteringinput-output data to remove the high frequency dither in the indoor airtemperature T_(ia) and can compute a time-averaged version of the HVACload {dot over (Q)}_(HVAC) from the discrete HVAC staging trajectory.The thermal model parameters and HVAC load model parameters θ can thenbe fit to the input-output data (or filtered input-output data). In someembodiments, step 1402 includes augmenting the resulting state-spacemodel with another integrating disturbance model and estimating theKalman filter gain for the resulting model. Using data obtained in asecondary experiment or under normal operation, step 1402 can includevalidating the model through the use of statistics that capture themulti-step prediction accuracy of the resulting model.

Still referring to FIG. 14, process 1400 is shown to include predictinga heat load disturbance {dot over (Q)}_(other) experienced by a buildingzone at each time step (k=0 . . . N−1) of an optimization period (step1404). In some embodiments, step 1404 is performed by load/ratepredictor 518 and state/disturbance predictor 520. Step 1404 can includegenerating a disturbance estimate {circumflex over (d)}(k) for each timestep k in the optimization period. In some embodiments, step 1404includes estimating the disturbance state {circumflex over (d)}(k) ateach time step k in the optimization period using the followingstate/disturbance model:

$\begin{bmatrix}{{\hat{T}}_{ia}\left( {k + 1} \middle| k \right)} \\{{\hat{T}}_{m}\left( {k + 1} \middle| k \right)} \\{\hat{I}\left( {k + 1} \middle| k \right)} \\{\hat{d}\left( {k + 1} \middle| k \right)}\end{bmatrix} = {{{A_{aug}\begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\hat{T}}_{m}\left( k \middle| {k - 1} \right)} \\{\hat{I}\left( k \middle| {k - 1} \right)} \\{\hat{d}\left( k \middle| {k - 1} \right)}\end{bmatrix}} + {B_{aug}\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}} + {\begin{bmatrix}K_{x} \\K_{d}\end{bmatrix}{\left( {\begin{bmatrix}{T_{ia}(k)} \\{{\overset{.}{Q}}_{HVAC}(k)}\end{bmatrix} - \begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\overset{.}{\hat{Q}}}_{HVAC}\left( k \middle| {k - 1} \right)}\end{bmatrix}} \right)\mspace{20mu}\begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\overset{.}{\hat{Q}}}_{HVAC}\left( k \middle| {k - 1} \right)}\end{bmatrix}}}} = {{C_{aug}\begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\hat{T}}_{m}\left( k \middle| {k - 1} \right)} \\{\hat{I}\left( k \middle| {k - 1} \right)} \\{\hat{d}\left( k \middle| {k - 1} \right)}\end{bmatrix}} + {D\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}}}}$$\mspace{20mu} {{{{where}\mspace{14mu} A_{aug}} = \begin{bmatrix}A & B_{d} \\0 & I\end{bmatrix}},{B_{aug} = \begin{bmatrix}B \\0\end{bmatrix}},{C_{aug} = {\begin{bmatrix}C & 0\end{bmatrix}.}}}$

Step 1404 can include recording a history of disturbance estimates{circumflex over (d)}_(hist):={{circumflex over (d)}(k|k−1)}_(k=1) ^(N)^(hist) predicted by the state/disturbance model.

Step 1404 can include using the history of disturbance estimates{circumflex over (d)}_(hist) to predict the value of {circumflex over({dot over (Q)})}_(other)(k) at each time step k. In some embodiments,the heat load disturbance {circumflex over ({dot over (Q)})}_(other)(k)is a function of the time of day, the day type (e.g., weekend, weekday,holiday), the weather forecasts, building occupancy, and/or otherfactors which can be used by load/rate predictor 518 to predict the heatload disturbance {circumflex over ({dot over (Q)})}_(other)(k). Forexample, step 1404 can include predicting the heat load disturbanceusing the equation:

{circumflex over ({dot over (Q)})}_(other)(k)=f _(pred)(kΔ+t ₀,{circumflex over (d)} _(hist))

where Δ>0 is the sample period and t₀ is an initial time. In someembodiments, step 1404 includes using weather forecasts from a weatherservice 526 and/or load history data from historical data 528 to predictthe heat load disturbance {circumflex over ({dot over (Q)})}_(other)(k)at each time step.

In some embodiments, step 1404 includes using a deterministic plusstochastic model trained from historical load data to predict the heatload disturbance {circumflex over ({dot over (Q)})}_(other)(k). Step1404 can include using any of a variety of prediction methods to predictthe heat load disturbance {circumflex over ({dot over (Q)})}_(other)(k)(e.g., linear regression for the deterministic portion and an AR modelfor the stochastic portion). In some embodiments, step 1404 includespredicting the heat load disturbance {circumflex over ({dot over(Q)})}_(other)(k) using the techniques described in U.S. patentapplication Ser. No. 14/717,593.

Still referring to FIG. 14, process 1400 is shown to include generatingan economic cost function that accounts for a cost of resources consumedby HVAC equipment at each time step (k=0 . . . N−1) of an optimizationperiod (step 1406). An example of an objective function which can begenerated in step 1406 is shown in the following equation:

$\min\limits_{T_{sp},\epsilon_{T},{\delta \; T_{sp}}}{\sum\limits_{k = 0}^{N - 1}\left( {{{r_{elec}(k)}{{\hat{}}_{elec}(k)}} + {p_{\epsilon_{T}}{\epsilon_{T}(k)}} + {p_{\delta \; T}\delta \; {T_{sp}(k)}}} \right)}$

where r_(elec)(k) is the price of electricity at time step k,

_(elec)(k) is the predicted electricity consumption of HVAC equipment308 at time step k, ∈_(T)(k) is a number of degrees by which the zonetemperature constraints are violated at time step k, p_(∈) _(T) is azone air temperature penalty coefficient applied to ∈_(T)(k), δT_(sp)(k)is a number of degrees by which the zone temperature setpoint T_(sp)changes between time step k−1 and time step k, and p_(δT) _(sp) is atemperature setpoint change penalty coefficient applied to δT_(sp)(k).Equation 24 accounts for the cost of electricity consumption(r_(elec)(k)

_(elec)(k)), penalizes violations of the indoor air temperature bounds(p_(∈) _(T) ∈_(T)(k)), and penalizes changes in the temperature setpoint(p_(δT)δT_(sp)(k)).

If additional resources (other than electricity) are consumed by HVACequipment 308, additional terms can be added to the objective functionto represent the cost of each resource consumed by HVAC equipment 308.For example, if HVAC equipment 308 consume natural gas and water inaddition to electricity, the objective function can be updated toinclude the terms r_(gas)(k)

_(gas)(k) and r_(water)(k)

_(water)(k) in the summation, where r_(gas)(k) is the cost per unit ofnatural gas at time step k,

_(gas)(k) is the predicted natural gas consumption of HVAC equipment 308at time step k, r_(water)(k) is the cost per unit of water at time stepk, and

_(water)(k) is the predicted water consumption of HVAC equipment 308 attime step k.

In some embodiments, step 1406 includes estimating the amount of eachresource consumed by HVAC equipment 308 (e.g., electricity, natural gas,water, etc.) as a function of the sensible heating or cooling load{circumflex over ({dot over (Q)})}_(HVAC). In some embodiments, aconstant efficiency model is used to compute the resource consumption ofHVAC equipment 308 as a fixed multiple heating or cooling load{circumflex over ({dot over (Q)})}_(HVAC) (e.g.,

_(elec)(k)=η{circumflex over ({dot over (Q)})}_(HVAC)(k)). In otherembodiments, equipment models describing a relationship between resourceconsumption and heating/cooling production can be used to approximatethe resource consumption of the HVAC equipment 308. The equipment modelsmay account for variations in equipment efficiency as a function of loadand/or other variables such as outdoor weather conditions (e.g.,

_(elec)(k)=f({circumflex over ({dot over (Q)})}_(HVAC)(k), T_(oa)(k),η_(base))). Step 1406 can include converting the estimatedheating/cooling load {circumflex over ({dot over (Q)})}_(HVAC)(k) ateach time step to one or more resource consumption values (e.g.,

_(elec)(k),

_(gas)(k),

_(water)(k), etc.) for inclusion in the objective function.

Still referring to FIG. 14, process 1400 is shown to include estimatinginitial system states {circumflex over (x)}(0) for a first time step(k=0) of the optimization period (step 1408). The initial system states{circumflex over (x)}(0) may include the estimated zone air temperature{circumflex over (T)}_(ia)(0), the estimated zone mass temperature{circumflex over (T)}_(m)(0), and the estimated integrating disturbanceÎ(0) at the beginning of the optimization period. In some embodiments,step 1408 is performed by state/disturbance estimator 520. The initialsystem states {circumflex over (x)}(0) can be estimated using same thestate/disturbance estimation model used in step 1404. For example, step1408 can include estimating the system states {circumflex over (x)}(k)at each time step k using the following equation:

{circumflex over (x)}(k+1|k)=A{circumflex over(x)}(k|k−1)+Bu(k)+K(y(k)−ŷ(k|k−1))

ŷ(k|k−1)=C{circumflex over (x)}(k|k−1)+Du(k)

where {circumflex over (x)}(k+1|k) is the estimated/predicted state attime step k+1 given the measurement at time step k, ŷ(k|k−1) is thepredicted output at time step k given the measurement at time step k−1,and K is the estimator gain.

In some embodiments, the estimated state vector {circumflex over(x)}(k+1|k), the output vector ŷ(k|k−1), and the input vector u(k) aredefined as follows:

${{\hat{x}\left( {k + 1} \middle| k \right)} = \begin{bmatrix}{{\hat{T}}_{ia}\left( {k + 1} \middle| k \right)} \\{{\hat{T}}_{m}\left( {k + 1} \middle| k \right)} \\{\hat{I}\left( {k + 1} \middle| k \right)}\end{bmatrix}},{{\hat{y}\left( k \middle| {k - 1} \right)} = \begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\overset{.}{\hat{Q}}}_{HVAC}\left( k \middle| {k - 1} \right)}\end{bmatrix}},{{u(k)} = \begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}}$

where {circumflex over (T)}_(ia) is an estimate of the zone airtemperature T_(ia), {circumflex over (T)}_(m) is an estimate of the zonemass temperature T_(m), Î is an estimate of the integrating disturbanceI, {circumflex over ({dot over (Q)})}_(HVAC) is an estimate of theheating or cooling load provided by HVAC equipment 308, T_(sp) is thetemperature setpoint, and T_(oa) is the outdoor air temperature. Theestimated system states {circumflex over (x)}(0) at the first time stepk=0 can be used as the initial system states.

Still referring to FIG. 14, process 1400 is shown to include using athermal mass storage predictive model to constrain system states{circumflex over (x)}(k+1|k) at time step k+1 to the system states (k)and the temperature setpoint T_(sp)(k) at time step k (step 1410). Anexample of a constraint which can be generated based on the thermal massstorage predictive model is:

$\begin{bmatrix}{{\hat{T}}_{ia}\left( {k + 1} \middle| k \right)} \\{{\hat{T}}_{m}\left( {k + 1} \middle| k \right)} \\{\hat{I}\left( {k + 1} \middle| k \right)}\end{bmatrix} = {{A\begin{bmatrix}{{\hat{T}}_{ia}(k)} \\{{\hat{T}}_{m}(k)} \\{\hat{I}(k)}\end{bmatrix}} + {B\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}} + {B_{d}{{\hat{\overset{.}{Q}}}_{other}(k)}}}$

where {circumflex over (T)}_(ia)(k+1) is the predicted temperature ofthe building zone at time step k+1, {circumflex over (T)}_(m)(k+1) isthe predicted temperature of the building mass at time step k+1, Î(k+1)is the predicted value of the integrating disturbance at time step k+1,{circumflex over (T)}_(ia)(k) is the predicted temperature of thebuilding zone at time step k, {circumflex over (T)}_(m)(k) is thepredicted temperature of the building mass at time step k, Î(k) is thepredicted value of the integrating disturbance at time step k, T_(sp)(k)is the temperature setpoint at time step k, T_(oa)(k) is the outdoor airtemperature (measured) at time step k, and {dot over ({circumflex over(Q)})}_(other)(k) is the estimated heat load disturbance at time step k.

Process 1400 is shown to include using a HVAC load predictive model toconstrain the HVAC equipment load {dot over ({circumflex over(Q)})}_(HVAC)(k) at time step k to systems states {circumflex over(x)}(k) and the temperature setpoint T_(sp)(k) at time step k (step1412). An example of a constraint which can be generated based on theHVAC load model is:

$\left\lbrack {{\hat{\overset{.}{Q}}}_{HVAC}(k)} \right\rbrack = {{C\begin{bmatrix}{{\hat{T}}_{ia}(k)} \\{{\hat{T}}_{m}(k)} \\{\hat{I}(k)}\end{bmatrix}} + {D\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}}}$

where {dot over ({circumflex over (Q)})}_(HVAC)(k) is the predicted HVACequipment load at time step k and {circumflex over (T)}_(ia)(k), Î(k),T_(sp)(k), and T_(oa)(k) are the same as the corresponding variables inthe thermal mass storage model.

Still referring to FIG. 14, process 1400 is shown to include optimizingthe economic cost function subject to constraints defined by thepredictive models and constraints on the predicted zone temperature{circumflex over (T)}_(ia)(k) to determine optimal temperature setpointsT_(sp) for the building zone for each time step (k=0 . . . N−1) of theoptimization period (step 1414). The objective function can be optimizedsubject to a set of constraints (e.g., Equations 25-34). The constraintsmay include the thermal mass storage model and the HVAC load model,which define the relationship between {circumflex over ({dot over(Q)})}_(HVAC), the temperature setpoints T_(sp), and the zone airtemperature T_(ia) at each time step. The constraints may also includeconstraints on the zone air temperature T_(ia)(k) and constraints thatdefine the penalty terms in the objective function. The optimizationperformed in step 1414 may generate a temperature setpoint T_(sp)(k) foreach time step in the optimization period.

Process 1400 is shown to include providing the temperature setpointT_(sp) for the first time step (k=0) of the optimization period to anequipment controller (step 1416). Step 1416 can include sending theoptimal temperature setpoints T_(sp) from model predictive controller302 to an equipment controller 406. In some embodiments, both modelpredictive controller 302 and equipment controller 406 are components ofa smart thermostat 100 (as shown in FIG. 4). In other embodiments, modelpredictive controller 302 is separate from smart thermostat 100 andconfigured to send the optimal temperature setpoints to smart thermostat100 via a communications network 304 (as shown in FIG. 3). In someembodiments, only the temperature setpoint for the first time step(i.e., T_(sp)(0)) is provided to the equipment controller in step 1416.The remaining temperature setpoints are for future time steps and may beupdated the next time the optimization is performed (e.g., at thebeginning of the next time step).

In some embodiments, step 1416 includes providing the entire set oftemperature setpoints T_(sp) (i.e., a temperature setpoint T_(sp)(k) foreach time step k) to the equipment controller. In the event thatcommunication between model predictive controller 302 and equipmentcontroller 406 is lost (e.g., network connectivity is disrupted),equipment controller 406 can use the set of temperature setpointsprovided in step 1416 until communications between model predictivecontroller 302 and equipment controller 406 are restored. For example,equipment controller 406 can use each temperature setpoint T_(sp)(k)received in step 1416 to control HVAC equipment 308 during thecorresponding time step k until communications between model predictivecontroller 302 and equipment controller 406 are restored.

Still referring to FIG. 14, process 1400 is shown to include usingfeedback collected from the HVAC system during the first time step (k=0)to update the thermal mass storage predictive model (step 1418). Thefeedback collected from the HVAC system may include measurements of thezone air temperature T_(ia) and/or measurements of the HVAC equipmentload {dot over (Q)}_(HVAC). Step 1418 may include calculating thedifference between the value of the zone air temperature {circumflexover (T)}_(ia)(k|k−1) predicted in step 1408 and the actual measuredvalue of the zone air temperature T_(ia)(k) collected in step 1418.Similarly, step 1418 may include calculating the difference between thevalue of the HVAC equipment load {circumflex over ({dot over(Q)})}_(HVAC)(k|k−1) predicted in step 1408 and the actual measuredvalue of the HVAC equipment load {dot over (Q)}_(HVAC)(k) collected instep 1418.

The updating performed in step 1418 may include multiplying the Kalmangain matrix K generated in step 1402 by the vector of differences, asshown in the following equation:

$\begin{bmatrix}{{\hat{T}}_{ia}\left( {k + 1} \middle| k \right)} \\{{\hat{T}}_{m}\left( {k + 1} \middle| k \right)} \\{\hat{I}\left( {k + 1} \middle| k \right)} \\{\hat{d}\left( {k + 1} \middle| k \right)}\end{bmatrix} = {{A_{aug}\begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\hat{T}}_{m}\left( k \middle| {k - 1} \right)} \\{\hat{I}\left( k \middle| {k - 1} \right)} \\{\hat{d}\left( k \middle| {k - 1} \right)}\end{bmatrix}} + {B_{aug}\begin{bmatrix}{T_{sp}(k)} \\{T_{oa}(k)}\end{bmatrix}} + {\begin{bmatrix}K_{x} \\K_{d}\end{bmatrix}\left( {\begin{bmatrix}{T_{ia}(k)} \\{{\overset{.}{Q}}_{HVAC}(k)}\end{bmatrix} - \begin{bmatrix}{{\hat{T}}_{ia}\left( k \middle| {k - 1} \right)} \\{{\overset{.}{\hat{Q}}}_{HVAC}\left( k \middle| {k - 1} \right)}\end{bmatrix}} \right)}}$

In some embodiments, the updating performed in step 1418 is equivalentto the updating step of a Kalman filter in which the predicted states{circumflex over (x)}(k+1|k) are calculated based on the error betweenthe predicted system output ŷ(k|k−1) at time step k and theactual/measured values of the predicted output variables y(k).

Process 1400 is shown to include shifting the optimization forward intime by one time step (step 1420) and returning to step 1408. Steps1408-1420 can be repeated at the beginning of each time step to generatea set of temperature setpoints T_(sp) for each time step in the shiftedoptimization period. The first temperature setpoint T_(sp)(0) generatedfor the shifted optimization period can be provided to the equipmentcontroller and used to control the HVAC equipment.

Performance Graphs

Referring now to FIGS. 15-18, several graphs 1500-1800 illustrating theperformance of model predictive control (MPC) systems 300-400 are shown,according to some embodiments. Graph 1500 compares the energyconsumption of MPC systems 300-400 to the energy consumption of abaseline temperature control system without MPC. Graph 1600 indicatesthe cost of electricity (i.e., $/kWh) which varies over time. Line 1502represents the energy consumption of the baseline system, whereas line1504 represents the energy consumption of MPC systems 300-400. MPCsystems 300-400 use more energy than the baseline system during anoff-peak period 1506 when energy prices are low. However, MPC systems300-400 use less energy than the baseline system during a peak period1508 when energy prices are high. Advantageously, this results in a costsavings relative to the baseline system.

Graph 1700 compares the building mass temperature trajectories of thebaseline system and MPC systems 300-400. Line 1702 represents thetemperature T_(m) of the solid mass within building zone 310 when thebaseline system is used to provide temperature control. Since thebaseline system does not store energy in the building mass, thetemperature T_(m) remains relatively constant over the duration of theoptimization period. Line 1704 represents the temperature T_(m) of thesolid mass within building zone 310 when MPC systems 300-400 are used toprovide temperature control. During the off-peak period 1506, buildingzone 310 is precooled, which results in a decrease in the temperatureT_(m) of the building mass. During the peak period 1508, thermal energyfrom the air within building zone 310 is moved into the building mass,which increases the temperature T_(m) of the building mass and providescooling for the air within building zone 310.

Graph 1800 compares the zone air temperature trajectories of thebaseline system and MPC systems 300-400. Line 1802 represents thetemperature T_(ia) of the air within building zone 310 when the baselinesystem is used to provide temperature control. The temperature T_(ia)remains relatively constant over the duration of the optimizationperiod. Line 1704 represents the zone air temperature T_(ia) withinbuilding zone 310 when MPC systems 300-400 are used to providetemperature control. During the off-peak period 1506, building zone 310is precooled to the minimum comfortable zone temperature T_(min) (e.g.,70° F.). The zone air temperature T_(ia) remains relatively constantduring the off-peak period 1506, but the temperature of the buildingmass T_(m) decreases as heat is removed from building zone 310. Duringthe peak period 1508, the building zone temperature T_(ia) is allowed toincrease to the maximum comfortable zone temperature T_(max) (e.g., 76°F.). Heat from the air within building zone 310 flows into the buildingmass, which provides cooling for the air within building zone 310 andincreases the temperature of the building mass T_(m).

Configuration of Exemplary Embodiments

The construction and arrangement of the systems and methods as shown inthe various exemplary embodiments are illustrative only. Although only afew embodiments have been described in detail in this disclosure, manymodifications are possible (e.g., variations in sizes, dimensions,structures, shapes and proportions of the various elements, values ofparameters, mounting arrangements, use of materials, colors,orientations, etc.). For example, the position of elements can bereversed or otherwise varied and the nature or number of discreteelements or positions can be altered or varied. Accordingly, all suchmodifications are intended to be included within the scope of thepresent disclosure. The order or sequence of any process or method stepscan be varied or re-sequenced according to alternative embodiments.Other substitutions, modifications, changes, and omissions can be madein the design, operating conditions and arrangement of the exemplaryembodiments without departing from the scope of the present disclosure.

The present disclosure contemplates methods, systems and programproducts on any machine-readable media for accomplishing variousoperations. The embodiments of the present disclosure can be implementedusing existing computer processors, or by a special purpose computerprocessor for an appropriate system, incorporated for this or anotherpurpose, or by a hardwired system. Embodiments within the scope of thepresent disclosure include program products comprising machine-readablemedia for carrying or having machine-executable instructions or datastructures stored thereon. Such machine-readable media can be anyavailable media that can be accessed by a general purpose or specialpurpose computer or other machine with a processor. By way of example,such machine-readable media can comprise RAM, ROM, EPROM, EEPROM, CD-ROMor other optical disk storage, magnetic disk storage or other magneticstorage devices, or any other medium which can be used to carry or storedesired program code in the form of machine-executable instructions ordata structures and which can be accessed by a general purpose orspecial purpose computer or other machine with a processor. Combinationsof the above are also included within the scope of machine-readablemedia. Machine-executable instructions include, for example,instructions and data which cause a general purpose computer, specialpurpose computer, or special purpose processing machines to perform acertain function or group of functions.

Although the figures show a specific order of method steps, the order ofthe steps may differ from what is depicted. Also two or more steps canbe performed concurrently or with partial concurrence. Such variationwill depend on the software and hardware systems chosen and on designerchoice. All such variations are within the scope of the disclosure.Likewise, software implementations could be accomplished with standardprogramming techniques with rule based logic and other logic toaccomplish the various connection steps, processing steps, comparisonsteps and decision steps.

What is claimed is:
 1. A system for monitoring and controllingtemperature of a building zone, the system comprising: a modelpredictive controller configured to: obtain a cost function thataccounts for a cost of operating HVAC equipment during each of aplurality of time steps in a period to provide heating or cooling to thebuilding zone; use a predictive model to predict the temperature of thebuilding zone during each of the plurality of time steps as a functionof temperature setpoints for a thermostat for the building zone; andgenerate the temperature setpoints for the thermostat for the buildingzone for each of the plurality of time steps by performing an operationof the cost function subject to a constraint on the temperature of thebuilding zone predicted by the predictive model; an equipment controllerconfigured to receive the temperature setpoints generated by the modelpredictive controller and drive the temperature of the building zonetoward the temperature setpoints during each of the plurality of timesteps by operating the HVAC equipment to provide heating or cooling tothe building zone.
 2. The system of claim 1, wherein the operation is anoptimization and the model predictive controller is configured todetermine the cost of operating the HVAC equipment during each of theplurality of time steps using a set of time-varying utility ratescomprising a utility rate value for each time step; wherein thetime-varying utility rates are received from a utility provider orpredicted by the model predictive controller.
 3. The system of claim 1,wherein the model predictive controller is separate from the thermostatand configured to provide the temperature setpoints to the thermostatvia a communications network.
 4. The system of claim 1, wherein theequipment controller is configured to: receive user-provided temperaturesetpoints via a local user interface of the thermostat or a mobiledevice in communication with the thermostat; and in response toreceiving the user-provided temperature setpoints, override thetemperature setpoints generated by the model predictive controller withthe user-provided temperature setpoints and use the user-providedtemperature setpoints to operate the HVAC equipment.
 5. The system ofclaim 1, wherein the model predictive controller is configured togenerate the predictive model by performing a system identificationprocess comprising: modulating the temperature setpoints within aconstrained temperature setpoint range; collecting a set of input-outputdata comprising values of the temperature setpoints and values of thetemperature of the building zone that result from modulating thetemperature setpoints during each of a plurality of time steps of alearning period; and fitting parameters of the predictive model to theset of input-output data.
 6. The system of claim 1, wherein thepredictive model comprises a thermal mass storage model that defines thetemperature of the building zone as a function of at least one of: heattransfer between air within the building zone and solid mass within thebuilding zone; heat transfer between the building zone and the HVACequipment; and an unmeasured heat load disturbance.
 7. The system ofclaim 6, wherein the model predictive controller is configured topredict a value of the unmeasured heat load disturbance experienced bythe building zone at each of the plurality of time steps in the period.8. The system of claim 1, wherein the predictive model comprises an HVACload model that defines the heating or cooling provided by the HVACequipment as a function of the temperature of the building zone and atleast one of the temperature setpoints.
 9. The system of claim 1,wherein the model predictive controller is configured to predict thecost of operating the HVAC equipment as a function of an amount of theheating or cooling provided by the HVAC equipment at each time step ofthe period.
 10. The system of claim 1, wherein the constraint on thetemperature of the building zone predicted by the predictive modelrequires the model predictive controller to maintain the temperature ofthe building zone predicted by the predictive model within: a first zonetemperature range during a first time step of the period; and a secondzone temperature range, different from the first zone temperature range,during another time step of the period subsequent to the first timestep.
 11. A method performed at least partially by a thermostat for abuilding zone for monitoring and controlling temperature of the buildingzone, the method comprising: obtaining a cost function that accounts fora cost of operating HVAC equipment during each of a plurality of timesteps in a period to provide heating or cooling to the building zone;using a predictive model to predict the temperature of the building zoneduring each of the plurality of time steps as a function of temperaturesetpoints for a thermostat for the building zone; generating thetemperature setpoints for the thermostat for the building zone for eachof the plurality of time steps by performing an operation of the costfunction subject to a constraint on the temperature of the building zonepredicted by the predictive model; operating the HVAC equipment toprovide heating or cooling to the building zone during each of theplurality of time steps to drive the temperature of the building zonetoward the temperature setpoints.
 12. The method of claim 11, whereinthe operation is an optimization, the method further comprising:receiving a set of time-varying utility rates from a utility provider orpredicting the time-varying utility rates, the set of time-varyingutility rates comprising a utility rate value for each time step; anddetermining the cost of operating the HVAC equipment during each of theplurality of time steps using the set of time-varying utility rates. 13.The method of claim 11, wherein the temperature setpoints are generatedby a model predictive controller separate from the thermostat andprovided to the thermostat via a communications network.
 14. The methodof claim 11, further comprising: receiving user-provided temperaturesetpoints via a local user interface of the thermostat or a mobiledevice in communication with the thermostat; and in response toreceiving the user-provided temperature setpoints, overriding thetemperature setpoints generated by the model predictive controller withthe user-provided temperature setpoints and using the user-providedtemperature setpoints to operate the HVAC equipment.
 15. The method ofclaim 11, further comprising generating the predictive model byperforming a system identification process comprising: modulating thetemperature setpoints within a constrained temperature setpoint range;collecting a set of input-output data comprising values of thetemperature setpoints and values of the temperature of the building zonethat result from modulating the temperature setpoints during each of aplurality of time steps of a learning period; and fitting parameters ofthe predictive model to the set of input-output data.
 16. The method ofclaim 11, wherein the predictive model comprises a thermal mass storagemodel that defines the temperature of the building zone as a function ofat least one of: heat transfer between air within the building zone andsolid mass within the building zone; heat transfer between the buildingzone and the HVAC equipment; and an unmeasured heat load disturbance.17. The method of claim 16, further comprising predicting a value of theunmeasured heat load disturbance experienced by the building zone ateach of the plurality of time steps in the period.
 18. The method ofclaim 11, wherein the predictive model comprises an HVAC load model thatdefines the heating or cooling provided by the HVAC equipment as afunction of the temperature of the building zone and at least one of thetemperature setpoints.
 19. The method of claim 11, further comprisingpredicting the cost of operating the HVAC equipment as a function of anamount of the heating or cooling provided by the HVAC equipment at eachtime step of the period.
 20. A system for monitoring and controllingtemperature of a building zone, the system comprising: a modelpredictive controller configured to: obtain a cost function thataccounts for a cost of operating HVAC equipment during each of aplurality of time steps in an optimization period to provide heating orcooling to the building zone; and generate temperature setpoints for thebuilding zone for each of the plurality of time steps by performing anoptimization of the cost function subject to a constraint on a predictedtemperature of the building zone, wherein the predicted temperature isdefined by a predictive model as a function of the temperaturesetpoints; a thermostat for the building zone, the thermostat configuredto: receive the temperature setpoints generated by the model predictivecontroller via a communications network; and drive the temperature ofthe building zone toward the temperature setpoints during each of theplurality of time steps by operating the HVAC equipment to provideheating or cooling to the building zone.